      Program MAIN_NEW
	USE cylinderGrid
	use accur
	IMPLICIT REAL(8)(A-H,K-M,O-Z)

      COMMON /BEGDATA/ rho0T,oblatnessT
      DOUBLE PRECISION rho0T,oblatnessT
      DOUBLE PRECISION J_CGS,J0,JK,JL,JR
      COMMON /GENRES1/ T1,T2,T3,M_S,J_CGS,TAUT,gammaT
     *       ,T8,T9,T10,T11,CM,T13,CJ,CI,T16,VT,T18
      COMMON /DOP/IOBHOD,IBT
      INTEGER INDEXT,IOBHOD,IB1,ILL,IRR,IND1
      Character*6 rotlaw

       open(1,file='rho0yudinm80.out',status='unknown')
       open(77,file='FLFR.out')
       write(1,'(18a)')' INDEX       rho0          ',
     *	'r(IB)**delta','     Re
     *	         ','       CM_u/S_mass ','     CJ_u/1d50 ','
     * T/abs(W) ','       gamma_  ',	'       rhomax  ','      phi0   '
     *,'        phiA   ','        phiB    ','        CM     ','
     *V     ','        CJ     ','        CI    ','         W     ',
     *'         VT     ','      Cpsi        IB'

C !!!  <<<<<<<<<<<<<<<<<<<<<<< 20.05.2018 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      dr=one/IA
      do i=0,IA
      r(i)=dr*i
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      enddo


!     rho0 is placed in module accur. SB: actually, in module cylinderGrid
!     Re is placed in cylinderGrid.

C --- ARRAY TTT CONTENTS THE NEXT INFORMATION:
!     TTT(1)=rho0T                                 =rho0
!	TTT(2)=r(IB)**delta                          =oblatnessT
!	TTT(3)=Re                                    =Re
!	TTT(4)=CM*Re**3*rho0/S_mass                  =M_S
!    	TTT(5)=CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50   =J_CGS
!     TTT(6)=T/abs(W)                              =TAUT
!	TTT(7)=Pip/Eth+one                           =gammaT
!	TTT(8)=rhomax
!	TTT(9)=phi(I0,JMax)
!     TTT(10)=phi(IA,JMax)
!	TTT(11)=phi(IB,JB)
!	TTT(12)=CM
!	TTT(13)=V
!	TTT(14)=CJ
!	TTT(15)=CI
!	TTT(16)=W
!	TTT(17)=VT                                   =VT
!	TTT(18)=Cpsi

!      IOBHOD=0
C ---- INDEXT - THROUGH ("SKVOZNOY") NUMBER OF COMPUTATION.
       INDEXT=1
	 rotlaw='rigid'
!      rotlaw='collap'
!      rotlaw='alpha'
!       rotlaw='omega'
C      rotlaw='saved'
! --- TWO NON STANDARD STEPS.
! --- THE FIRST NON STANDARD STEP. ----------------------
!    		rho0T=3.00D14
!      	Drho0=0.0050D14
!     Drho0R=0.00125*10.  ! RELATIVE VARIATION OF DENSITY IN CENTRE.
!       rho0T=3.00015752E+14 ! 0.95 continue
!       rho0T=3.12113244E+14 ! 0.9 continue
!      rho0T=3.31182278E+14 ! 0.9 - continue 0.7
CC        rho0T=3.32846510E+14 ! Begin M=0.2
      rho0T=3.D+14
!       rho0T=2.24914063E+14 ! Begin M=0.11

       Drho0R=0.005D0  ! It was 0.01,0.005.0025,0.00125. 23.01.2019.
       Drho0=rho0T*Drho0R

!   	rho0T=rho0T*(1.D0-Drho0R)

!      DDRFIX=0.05  ! 5% ! relative variation of Re.  25.08.2018.
!	 DDRFIX=0.1   ! 10%
	 DDRFIX=0.05  ! 1%, 2%  ! It was 0.07. 23.01.2019.
	 IOBHOD=0

!	OBLATNESS=0.9
      oblatnessT=0.7
      IBT=oblatnessT*IA+0.5
	oblatnessT=r(IBT)

      Print*,'IB-- !!!!==== IB,',IBT,IA,oblatnessT
!      READ*


	CALL rho0Yudinnmain(1,rotlaw)
      J0=J_CGS
      M0=M_S
      CM0=CM
      CI0=CI
      Print*,'First', J0,M0,CM0,CI0
      READ*
	STOP

	RR1=Re
      IND1=100
!      rho0T=rho0T-Drho0
!  	rho0T=rho0T*(1.D0-Drho0R)

	IBT=IBT+1

      OPEN(101,file='RESULT.dat')
	DO
	CALL rho0Yudinnmain(10,rotlaw)
      JK=J_CGS
      MK=M_S
      CMK=CM
      CIK=CI
	print*,rho0T,IBT
       Print*,'Next ',JK,MK,CM,CI
	 Print*,'Zero ',J0,M0,CM0,CI0
!		 	 Read*
      FL=-(JK-J0)/(0.5D0*(JK+J0))
   	FR=-(MK-M0)/(0.5D0*(MK+M0))*(CMK+CM0)/(CIK+CI0)
	FRR=-(MK-M0)/(MK+M0)*(CMK/CIK+CM0/CI0)
	print*,FL,FR
	WRITE(101,'(20ES23.15)'),rho0T,JK,J0,MK,M0,CM,
     *	CM0,CI,CI0,FL,FR,FRR
	read*
	rho0T=rho0T*1.00001D0
	END DO



      I2=1
!        *************************************************************
          DO I2=2,500   ! THIS IS EXTERNAL CYCLE ALONG THE DENSITY.
	 IOBHOD=1
!	 IOBHOD=0
 	CALL rho0Yudinnmain(I2,rotlaw)
      JK=J_CGS
      MK=M_S
      CMK=CM
      CIK=CI
       Print*,'**** MK,M0,CM,CI',MK,M0,CM,CI,rho0T
!		 	 Read*
      FL=-(JK-J0)/(0.5D0*(JK+J0))
   	FR=-(MK-M0)/(0.5D0*(MK+M0))*(CMK+CM0)/(CIK+CI0)

      WRITE (77,*) T2,FL,FR !,rho0T
!	WRITE (77,*) T2,FL,FR,rho0T,IBT,JK,J0,MK,M0,CMK,CM0,CIK,CI0
!	     Print*,'**** J0,JK,FL,FR,IBT',J0,JK,FL,FR,IBT
!		 	 Read*
!	IF (JK.GE.J0) GO TO 12
! ----- GENERAL BRANCH.
         DO I1=2,IBT-20  ! CYCLE ON OBLATNESS.
      IND1=IND1+1
 !     IOBHOD=1 !!!!! 1
!      IOBHOD=0
!      IBT=IBT-1 !!!
	 IBT=IBT+1

       IF (IBT.GE.(IA-1)) GO TO 11
!	J0=JK
!	M0=MK
!     CM0=CMK
!     CI0=CIK
	CALL rho0Yudinnmain(IND1,rotlaw)
      JK=J_CGS
	MK=M_S
      CMK=CM
      CIK=CI

	FL=-(JK-J0)/(0.5D0*(JK+J0))
	FR=-(MK-M0)/(0.5D0*(MK+M0))*(CMK+CM0)/(CIK+CI0)

       WRITE (77,*) T2,FL,FR  !,rho0T
!      WRITE (77,*) T2,FL,FR,rho0T,IBT,JK,J0,MK,M0,CMK,CM0,CIK,CI0
       Print*,'**** JK,J0,MK,M0,CM,CI',JK,J0,MK,M0,CM,CI,rho0T
	 Print*,'**** FL,FR,IBT',FL,FR,IBT
!		 	 Read*
!	IF (JK.GE.J0) GO TO 10 ! ZADANIE 1.
!	IF (FL.LE.FR) GO TO 10 ! ZADANIE 2. WITH UTOCHNENIE.
	IF (FL.GE.FR) GO TO 30 ! ZADANIE 2. WITHOUT UTOCHNENIE. 16.02.
         END DO
      GO TO 11
!	     Print*,'ERROR !!!!**** J0,IBT',J0,IBT
!		 	 Read*
! --- ------------------------
   10 CONTINUE

! ---- UTOCHNENIE ------09.01.2019.
      ITERM=20
      EPSUT=1.D-4
!      Drho0RUT=0.01  ! IT WAS.
       Drho0RUT=0.02
       IND2=200
       ITER=0
      rho0TUT= rho0T
! ---- BEGIN OF ITERATION
   21 CONTINUE
      ITER=ITER+1
	IF (ITER.LE.ITERM) GO TO 22
      Print*,'NUMBER OF ITERATIONS BECAME GREATER ITERM',rho0T,J0,IBT
		 	 Read*
      GO TO 11
   22 CONTINUE

      IND2=IND2+1
      A=rho0T
      CALL rho0Yudinnmain(IND2,rotlaw)
      JK=J_CGS
      FA=JK-J0

      IND2=IND2+1
	 rho0T=rho0T*(1.D0-Drho0RUT)
      B=rho0T
      CALL rho0Yudinnmain(IND2,rotlaw)
       JK=J_CGS
       FB=JK-J0
! ------- CHECKING OF SIGN.
      IF (FA*FB.LE.0.D0) GO TO 23
      Print*,'!!! LET US INCREASE Drho0RUT - A,B,FA,FB', A,B,FA,FB
		 	 Read*
      STOP
   23 CONTINUE
      Print*,'!!!!!!!! ---- A,B,FA,FB', A,B,FA,FB
!		 	 Read*

   24 CONTINUE
      XX=(A+B)/2.D0
      IND2=IND2+1
	 rho0T=XX
       CALL rho0Yudinnmain(IND2,rotlaw)
       JK=J_CGS
       FX=JK-J0

	IF ((FA*FX).LT.0.D0) THEN
	    B=XX
      ELSE
	   IF ((FB*FX).LT.0.D0) A=XX
!!!    ?????????
	END IF

! ----- CHECK OF FINISH OF ITERATIONS.
      IF ((DABS(A-B)/(0.5*(A+B))).GT.EPSUT) GO TO 24
      Print*,'FINISH ---- A,B,FA,FB', A,B,FA,FB
!		 	 Read*
! --- ------------------------------------------------

! ---   WRITING IN GENERAL MASSIV.
	 rho0T=XX

   30 CONTINUE ! WITHOUT  UTOCHNENIE 16.02.2019.
      IOBHOD=0
	CALL rho0Yudinnmain(I2,rotlaw)

      J0=J_CGS
      M0=M_S
      CM0=CM
      CI0=CI

! ---   PREPARING FOR NEW STEP ALONG DENSITY.
       DDRE=(DABS(Re-RR1))/0.5D0/(Re+RR1)
       IF (DDRE.GT.DDRFIX) Drho0R=Drho0R*0.7D0
          RR1=Re
!      	rho0T=rho0T-Drho0
  	rho0T=rho0T*(1.D0-Drho0R)
       END DO

      GO TO 11
! --- ------------------------
   12 CONTINUE
	     Print*,'ERROR * J0,JK,IBT',J0,JK,IBT ! make new print!
		 	 Read*
! --- ------------------------
! --- ------------------------
   11 CONTINUE

      CLOSE(77)
      CLOSE(1)
	Stop
	End
! ******************************************************************
! ******************************************************************

!         ______________________________________________________
!         ______________________________________________________
!         ______________________________________________________




C !!! 21.05.2018.
C- f for sm calls on HPUX and AIX, use '' for others

	SUBROUTINE rho0Yudinnmain(INDEXT,rotlaw)
C !!!       Program rho0Yudinnmain
C- _TRACE "@wterm' IB=',IB,' Ncheb=',Ncheb,' nDersMAX=',nDersMAX," ;
C- uses eosyudin.trf
C- constructs stationary barotropic configurations in fast rotation along rho0 series
C- this code may be used for testing non-rotating models
C- in this test one can put delta==1, while normally for alpha-law delta==2
C- but now delta=2 is not plotted correctly
C- the method is described in AA290(1994)674 for ursos
C- main variables:
C- phi(*,*) - grav.potential \hat\Phi as in eq.2 and 10
C- Fpsi - centrifugal potential \hat\Psi as in eq.8 and 14
C-*****************************************************************
      use utilities_module
C      use,intrinsic :: iso_c_binding
C      use dflib
      use cylinderGrid
      use accur
C-        _include accur;
C-        include 'accur.for';
C-        double precision  pi;
C:  for some functions *
      interface
      doubleprecision Function Fpsi(rcylw,rotlaw)
      implicit none
      doubleprecision::rcylw
      character*(*)::rotlaw
      endfunction
      doubleprecision Function Fdpsi_drcyl(rcylw,rotlaw)
      implicit none
      doubleprecision::rcylw
      character*(*)::rotlaw
      endfunction
      doubleprecision Function chi(rcylw,rotlaw)
C- \hat\Psi^0
      implicit none
      doubleprecision::rcylw
      character*(*)::rotlaw
      endfunction
      doubleprecision Function powrsmart(x,alpha)
      implicit none
      doubleprecision::alpha,x
      endfunction
      endinterface
C-       include 'statd.for';
C-        _include grid;
      doubleprecision S_mass,G_N,urho,uH
      Parameter(S_mass=1.989d33,G_N=6.6732d-8)
      Parameter(urho=2d9,uH=1.7d18)
C- WD
C: parameters like S_mass, Imax0 and commons *
C-       Parameter(S_mass=1.989d33, G_N=6.6732d-8); -- now in WDunits.inc
C       integer NrDAT,NzDAT,NtDAT,indmin(1),indmax(1);
C       Parameter(NrDAT=IA,NzDAT=4*IA,NtDAT=JMax);
C       REAL X(NrDAT+1),Y(NzDAT+1),phiRcylZplot(NrDAT,NzDAT);
C--    NrDAT=IA;   -- points on Rcylindrical
C--    NzDAT=4*IA; -- points on Z -- more points better for flattened bodies
C--    NtDAT=JMax; -- points on thet
C       integer NDIMP,NDP,NIP,ir,jt,i,j, -- for spherical coords
C               icyl,jcyl; -- for cylindr. coords
C       Parameter(irho0max=@irho0max);
C       double precision  rspl(0:NrDAT),thet(0:NtDAT),Phirthet(0:NrDAT,0:NtDAT),
C       rhoRcylZ(0:NzDAT,0:NrDAT)/,rhodV(0:NzDAT)*
C-     E01DAF  based on Example Program Text
C-     Mark 14 Release.  NAG Copyright 1989.
C-     .. Parameters ..
C-      INTEGER          NIN, *;
C-       PARAMETER        (NIN=5,*=6);
      INTEGER M,LIWRK,LWRK
      PARAMETER(LIWRK=NrDAT+1+2*(NrDAT-2)*(NtDAT-2),
     *	LWRK=(NrDAT+7)*(NtDAT+7))
      INTEGER ChMMAX,ChNMAX,nDerMX,ChLWRK,ChLIWRK
      PARAMETER(ChMMAX=NrDAT+1,ChNMAX=NrDAT+1,
C- ChNMAX==ChMMAX when no derivatives are given
     *nDerMX=2,ChLWRK=7*ChNMAX+5*nDerMX+ChMMAX+7,ChLIWRK=2*ChMMAX+2)
      INTEGER Ncheb,nDers(ChMMAX),skipCheb
      DOUBLEPRECISION Xch(ChMMAX),Ych(ChNMAX)
      DOUBLEPRECISION chebA(ChNMAX),ChWRK(ChLWRK),resCheb(ChMMAX),
     *	pt(ChMMAX)
C-     .. Local Scalars ..
C-       DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO;
      INTEGER IFAIL,MX,MY,PX,PY
C-     .. Local Arrays ..
      DOUBLEPRECISION Cspl((NrDAT+1)*(NtDAT+1)),Fspl((NrDAT+1)*(NtDAT+
     *	1)),FF(NzDAT+1),LAMDA(NrDAT+5),MU(NtDAT+5),WRK(LWRK),Xrad(0:
     *    NzDAT),Ythet(0:NzDAT)
      INTEGER IWRK(LIWRK)
C-     .. External Subroutines ..
      EXTERNAL E01DAF,E02DEF
C-     .. Intrinsic Functions ..
      INTRINSIC MAX,MIN
      Integer mmax,Imax0,ii,im,ip,NtableIB,nDersMAX,nDersMX,nDers1,Narg,
     *	KMax,nnzero,LICN,LIRN,NOUT,L,irho0,IY,N_max,nStep,Mcheb,ITMIN
     *    ,ITMAX,IRES,krho,itable
      doubleprecision Xmin,Xmax,Zcyl,CMGM,expand,x_min,x_max,ploty_min,
     *	ploty_min2,ploty_max,ploty_max2,plotx,ploty,ploty2
      parameter(mmax=1000)
      dimension plotx(0:mmax),ploty(0:mmax),ploty2(0:mmax)
C- for sm
      character*80 filein,psfile,device,titlestr,xlabel,ylabel,string,
     *	charg
      Parameter(IMax0=1000)
      doubleprecision EOS,rcylm,rcylp,xm,xp,C,C0,c1,CM,V,dV,CJ,CI,T,W,
     *	Eth,Pip,VT,rhomax,H,P,psi,psiOut,omegaA
C-        double precision r, dr, dr0, Re, alpha, dalpha, theta, dtheta, phi, phi0, phiu, phiA,
C-          epsMain, rho, vphi, vphi2, Fdpsi_drcyl;
C-        Dimension r(0:IMax),alpha(0:IMax),theta(0:JMax);
C-      real etime,dtime,tarray(2) ! Convex
C-        Common /phi/phi(0:IMax,0:JMax);
C-        double precision rho0,rho01,rho02,delta,gamma,CN,CK;
C-        Common /delta/delta;
C-        Common /gamma/gamma,CN,CK;
C-        Common /pi/pi;
C-        Common /IB/IB/JB/JB/I0/I0;
      doubleprecision phi0
      Dimension phi0(0:IMax0)
	doubleprecision myrotlaw,X1,X2
      Parameter(NtableIB=6)
      Real*8 tableIB(NtableIB)
C      data tableIB/1.,.917,.833,.750,.667,.625,.620,.619,.618,.617/
	data tableIB/0.34D0,0.339D0,0.338D0,0.337D0,0.336D0,0.335D0/
C      data tableIB/0.2D0,0.201D0,0.202D0,0.203D0,0.204D0,0.3D0,0.32D0,
C     *    0.34D0,0.37D0,0.4D0,0.45D0,0.5D0,0.6D0,0.7D0,0.8D0/
C- Hachisu
C-       Real*8 tableIB(NtableIB)/1.,.9,.8,.75,.7,.65,.64,.635,.63,.625/; -- AA290(1994)674
      Logical LogIter,LogInit,LogRapid,Lphi,Lrho,Lpolytrope
      Character*80 str,choblatness
      Character*6 rotlaw
C- rigid, jmscf, vconst, jconst or alpha
      REAL*4 secnds,t0,t1

! <<<<<<<<<<<<<<<<<<<<<<<  14.06.2018.
       COMMON /BEGDATA/ rho01T,oblatnessT
	 DOUBLE PRECISION rho01T,oblatnessT
	 DOUBLE PRECISION TTT
       COMMON /GENRES1/ TTT(1:18)
       INTEGER INDEXT,IOBHOD ! line 760
	 COMMON /DOP/IOBHOD,IBT
c >>>>>>>>>>>>>>>>>>>>>>>

C- timer
C-       LogRapid=.true.; -- default for CPsi in new_phi -- alpha-law rotation
      LogRapid=.false.
C- must be for other rotation laws test
C:  read rotlaw to be used  *
C------- '-->Entering Node %_arguments:'
      write(*,*)' Arguments needed: rotlaw -- either rigid,
     *	jmscf, vconst, jconst, alpha, collap, omega'
      write(*,*)' 2nd argument: oblatness: 0.< oblatness <= 1. '
      Narg=Iargc()
      GOTO(09999,09998,09997),Narg+1
      write(*,*)' Extra arguments.'
      Stop 16
      GOTO 09996
09999 CONTINUE ! Here the arguments are fixed
!      rotlaw='rigid'
C      rotlaw='collap'
C      rotlaw='alpha'
C      rotlaw='omega'
C      rotlaw='saved'
c !!!     oblatness=0.6d0  ! OBLATNESS!!!!!!!!!!!!!!!!!!!!!!!!
      oblatness=oblatnessT
      GOTO 09996
09998 CONTINUE
      call GetArg(1,rotlaw)
      oblatness=0.8d0
      GOTO 09996
09997 CONTINUE
      call GetArg(1,rotlaw)
      oblatness=1.d0
      call GetArg(2,choblatness)
      read(choblatness,*)oblatness
09996 CONTINUE

C !!!  <<<<<<<<<<<<<<<<<<<<<<< 20.05.2018 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!!!            dr=one/IA
!!!      do i=0,IA
!!!      r(i)=dr*i
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
!!!      enddo

C ---  A NEW VALUE OF oblatness in accodance with points (uzlami)
C      of grid along r.
       oblatness_old=oblatness
	 Print*,oblatness
      eps_ta=1.d-7
	 DO i=1,IA
         IF (((r(i-1)/r(IA)+eps_ta).LT.oblatness_old).AND.
     *	   ((r(i)/r(IA)+eps_ta).GE.oblatness_old)) then
            oblatness=r(i)/r(IA)
	   endif

	 Print *, oblatness,r(i), r(IA),i
	 ENDDO

	Print *, oblatness,oblatnessT,oblatnessA, r(IA),i
!	read*

C !!!  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      write(*,'(3a)')'  rotlaw=',rotlaw(1:len_trim(rotlaw)),' OK?'
      write(*,'(a,1p,1g11.3,a)')' oblatness=',oblatness,' OK?'
!      call mypause
      if(rotlaw(1:len_trim(rotlaw)).EQ.'alpha')then
      LogRapid=.true.
      else
      LogRapid=.false.
C- LogRapid is always true when rotlaw is not alpha
      endif
C------- '<--Leaving  Node %_arguments:'
C- control vars, better read them from a dat file:
      JB=0
C- spheroidal configurations
      Lrho=.true.
C- the set of configurations along rho0 sequence
ccc      Lrho=.false. ! -- the set of configurations along IB  sequence
      Lphi=.true.
C- initiate phi
C-      Lphi=.false.; -- read phi from file phi
C: constants and (r,theta) grid *
C------- '-->Entering Node %_init:'
      NOUT=16
c !!!      open(1,file='rho0yudinm80.out',status='unknown')
      open(2,file='f.dat')
      open(3,file='inputSeb.dat')
      open(4,file='rho0MJ.cur',access='append')
      open(7,file='rCJMccyl.out',status='unknown')
C- this is for rho centered
      open(8,file='rCJMcylgrid.out',status='unknown')
C- this is for rho on cyl.grid
      open(NOUT,file='rho0Cheb.out',status='unknown')
      open(12,file='Prseb.out',status='unknown')
C SB added
C- only for alpha-law
      open(13,file='alphaOmega.out',status='unknown')
      read(3,'(f15.10,a80)')delta
C- grid parameter 1 or 2, true radius is r**delta,
C- Aksenov used delta=2
      read(3,'(f15.10,a80)')gamma
C-=4d0/3d0;
C-       gamma=5d0/3d0;
C-        gamma=2d0;
      CN=one/(gamma-one)
      CK=one/(one+CN)
C-       rho0=2d14; -- for neutron stars
      read(3,'(e20.10,a80)')rho0
C-=2.d09; -- for polytropic and eosseb
      read(3,'(L20,a80)')Lpolytrope
C- true for polytropic case
      close(3)
      LinitMcyl=.true.
C- mcyl must be initiated for jmscf rotating law
      write(*,'(a,1p,3g15.8)')'delta,gamma,rho0',delta,gamma,rho0
C-        pause;
C-      Call timer(ITime0) ! PC
      t0=secnds(0.0)

      write(2,'(3a)')'         rho0,           r(IB)**delta,      Re,','
     *	             Mass,             V,           CJ,
     *  CI,','       T/abs(W),       Pip/Eth+1,         VT'
      write(12,'(7a)')'#    r(i)   ','  rho(i,0)  ','     P      ',
     *	'     H      ','  phi(i,0)  ','    psi     ','  vphi'

C !!!      dr=one/IA
C !!!      do i=0,IA
C !!!      r(i)=dr*i
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
C !!!      enddo
      dalpha=one/r(IA-1)-one/r(IA)
      do l=0,IMax-IA+1
      alpha(l)=dalpha*l
      enddo
      do i=IA+1,IMax-1
      l=IMax-i
      r(i)=one/alpha(l)
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      enddo
      dtheta=half*pi/dble(JMax)
      do j=0,JMax
      	theta(j)=dtheta*j
      enddo
      do ir=0,NrDAT
      Rcyld(ir)=powrsmart(r(ir),delta)
      enddo
C------- '<--Leaving  Node %_init:'
C:   initial phi - potential  *
C------- '-->Entering Node %_initphi:'
C-     === initial phi ==>
      if(Lphi)then
      do j=0,JMax
      do i=0,IMax
      if(i.le.IA)then
      phi(i,j)=float(i-IA-1)/float(IA)
      else
      l=IMax-i
      phi(i,j)=phi(IA,j)*powrsmart(alpha(l),delta)
      endif
      enddo
      enddo
      else
      open(16,file='phi',status='old')
      read(16,*)phi0
      close(16)
      do j=0,JMax
      do i=0,IMax
      if(i.le.IA)then
      dr0=one/IMax0
      im=powrsmart(r(i),delta)/dr0
      ip=im+1
      rcylm=im*dr0
      rcylp=ip*dr0
      xm=(rcylp-powrsmart(r(i),delta))/dr0
      xp=(powrsmart(r(i),delta)-rcylm)/dr0
      phi(i,j)=xm*phi0(im)+xp*phi0(ip)
      else
      l=IMax-i
      phi(i,j)=phi(IA,j)*powrsmart(alpha(l),delta)
      endif
      enddo
      enddo
      endif
C-     <== initial phi ===
C------- '<--Leaving  Node %_initphi:'
Cinitmcyl : initial mcyl()  *
      if(JB.eq.0)then
C- spheroidal configurations
      I0=0
      else
C- toroidal configurations
      I0=half*(IB+IA)
      endif
      if(Lrho)then
       rho01=rho01T
       rho02=rho01
	rho0=rho01   !!!
	Print*,' !!!!---- rho0,rho01', rho0,rho01
!	Read*

c      rho01=3.D14   !   DEN1----------------------------------------------------------------------
c      rho02=3.D14
c      rho02=1.5D14   !   DEN2----------------------------------------------------------------------
C- NS -- add to arguments
C-            rho01=5d4;
C           rho01=1d9;
C           rho02=5d9;  -- WD*
C-            irho0max=60; -- 20; -- not more than 999 for correct names of plot files
C-            write(4,'(a)')'  rho//J            M';
c !!!      do irho0=0,irho0max
       do irho0=0,0
C- Aksenov: ===== rho0 cycle ====>
C-                IB=IA; -- if IB==IA we have spherical configurations
C- IB=IA-80; -- if IB==IA we have spherical configurations
C- how to make non-zero rotation?
C-               IB=0.6*IA; -- try rotating sequence along rho0
!  <<<<<<<<<

      IB=oblatness*IA
      IB=IBT
	AIB=IB
	AIA=IA
      oblatness=AIB/AIA
      	print*,'OOOOO=',IBT,IB,AIB,oblatness,r(IB)
!	read*  !
!  >>>>>>>>>

!C- assuming irho0<=9999
!      IF(.NOT.(irho0.LT.10))GOTO 09995
!      write(str,'(a,I1)')'000',irho0
!      GOTO 09992
!09995 CONTINUE
!      IF(.NOT.(irho0.LT.100))GOTO 09994
!      write(str,'(a,I2)')'00',irho0
!      GOTO 09992
!09994 CONTINUE
!      IF(.NOT.(irho0.LT.1000))GOTO 09993
!      write(str,'(a,I3)')'0',irho0
!      GOTO 09992
!09993 CONTINUE
!      write(str,'(I4)')irho0
!09992 CONTINUE

!     <<<<<<<<<<<<<<<<<<<<<<<<< 17.06.2018.
C- assuming INDEXT<=9999
      IF(.NOT.(INDEXT.LT.10))GOTO 09995
      write(str,'(a,I1)')'000',INDEXT
      GOTO 09992
09995 CONTINUE
      IF(.NOT.(INDEXT.LT.100))GOTO 09994
      write(str,'(a,I2)')'00',INDEXT
      GOTO 09992
09994 CONTINUE
      IF(.NOT.(INDEXT.LT.1000))GOTO 09993
      write(str,'(a,I3)')'0',INDEXT
      GOTO 09992
09993 CONTINUE
      write(str,'(I4)')INDEXT
09992 CONTINUE

!     >>>>>>>>>>>>>>>>>>>>>>>>>
      open(26+INDEXT,file='phirho'//str(1:len_trim(str))//'.res')
      rho0=rho01 !!!+dble(irho0)/max(one,dble(irho0max))*(rho02-rho01)
      alphaPolytrope=sqrt((one+CN)*CK*rho0**(1.d0/CN-1.d0)/
     *	(4.d0*pi*G_N))
      write(*,'(a,1p,g12.3)')' rho0=',rho0
C-                call mypause;
C: Iterate Phi, C  *
C------- '-->Entering Node %_FindPhi:'
      write(*,*)' entering _FindPhi'
C-            call mypause;
      LogInit=.true.
      write(*,'(a,1p,g15.5,i7)')'rho0, IB',rho0,IB
      write(*,'(a,1p,3g15.5)')'phi(IA,JMax),phi(I0,JMax) =',phi(IA,JMax)
     *	,phi(I0,JMax)
C-           'phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)))=',phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)));
C- write of Fpsi is forbidden if there is write in it
      write(*,'(2a)')' rotlaw=',rotlaw(1:len_trim(rotlaw))
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)*
     *	EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      write(*,'(a,1p,3g15.5)')'phiu,EOS(rho0,3),C=',phiu,EOS(rho0,3),C
C-            call mypause;
      n_max=300
!     epsMain=1.d-12
!      epsMain=1.d-9
       epsMain=1.d-8
C- to data file
      LogIter=.true.
      nStep=1
      do while(LogIter.and.(nStep.LE.n_max))
      C0=C
      Call new_phi(
C-r,alpha,theta,dtheta, -- now in module cylinderGrid
     *C,LogInit,LogRapid,rotlaw(1:len_trim(rotlaw)))
C- last argument is the rotation law
C              if(C>C0)then;
C                C=C0+min(abs(C-C0),0.03d0*abs(C0));
C              else;
C                C=C0-min(abs(C-C0),0.03d0*abs(C0));
C              endif;*
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18

!       pRINT*, '========1'
!	READ*

      if(Dabs((C0-C)/C).LT.epsMain)LogIter=.false.
      write(*,*)' nStep,C,phi(I0,JMax):',nStep,C,phi(I0,JMax)
      nStep=nStep+1
      enddo
C------- '<--Leaving  Node %_FindPhi:'
C: to channel 26+irho0 after convergence, find integrals,
C                                remap onto cylindrical coordinates,
C                                prepare plots etc. *
C------- '-->Entering Node %_OutputResults:'
      Call Integrals(CM,V,CJ,CI,T,W,Eth,Pip,C,omegaA,
     *	rotlaw(1:len_trim(rotlaw)))
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' VT=',VT
      rhomax=zero
      do i=0,IA
C- loop on r
      do j=0,JMax
C- loop on theta
      if(j.EQ.JMax)then
      phiEqu(i)=phi(i,j)
C                   write(*,'(a,2i5,1p,2g15.6)')' i, j, theta(j),phiEqu(i)=',
C                            i,j,theta(j),phiEqu(i);
C                   if(mod(i,20)==0)call mypause;*
      endif
C- write(*,*)i,j,delta,C;
      rho=EOS(phiu*(C-phi(i,j)-Fpsi(powrsmart(r(i),delta)*sin(theta(j))
     *	,rotlaw(1:len_trim(rotlaw)))),1)
C- in phys. units from EOS
      rhoOn_r_thetaGrid(i,j)=rho
C- save for future
      rhomax=max(rhomax,rho/rho0)
      if(i.GT.0.and.j.EQ.Jmax.and.irho0.EQ.irho0max)then
      H=EOS(rho,3)
      P=EOS(H,4)
      psi=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
      vphi2=-Fdpsi_drcyl(powrsmart(r(i),delta)*sin(theta(j)),
     *rotlaw(1:len_trim(rotlaw)))*powrsmart(r(i),delta)*sin(theta(JMax))
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,vphi2))
      write(12,'(1p,7g12.4)')r(i),rho,P,H,phi(i,j),psi,vphi
      endif
C: write to channel 26+irho0 for plots and debug
C                              in spherical coordinates  *
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      psiOut=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
C-                  Phirthet(i,j) = C -phi(i,j)-psiOut;
      Phirthet(i,j)=phi(i,j)
      write(26+INDEXT,'(1p,3e20.12)')powrsmart(r(i),delta),theta(j),
     *	C-phi(i,j)-psiOut

CC SB added
      if(j.EQ.JMax)then
      if(i.EQ.0)then
      write(13,'(a,i5)')'# IB=',IB
      write(13,'(6(a,9x))')'#  i','  r','alpha1','alpha2','gradPhi',
     * 	'Omega^2'
      else
C                     write(*,*) 'i,r(i),theta(j)=',i,r(i),theta(j);
C                     write(*,'(a,2i5,1p4g12.3)')'i j phi Fpsi Fpsi/phi:',
C                            i,j,phi(i,j), psiOut, psiOut/(phi(i,j)-phi(I0,JMax));
C                            -- Phi0 == Phi_0 == phi(I0,JMax)
C                     write(*,'(a,1p4g12.3)')'phiB, phiA, Phi0, (phiB-phiA)/(PhiA-Phi0) :',
C                            phi(IB,0), phi(IA,JMax), phi(I0,JMax),
C                            (phi(IB,0)-phi(IA,JMax))/(phi(IA,JMax)-phi(I0,JMax));
C                     pause;
C                            *
      alpha1=psiOut/(phi(i,j)-phi(I0,JMax))
      alpha2=(phi(IB,0)-phi(IA,JMax))/(phi(IA,JMax)-phi(I0,JMax))
      gradPhi=(phi(i+1,Jmax)-phi(i-1,Jmax))/(2d0*dr)
      write(13,'(i5,1p,4g15.5,g17.5)')i,r(i),alpha1,alpha2,gradPhi,
     *	-alpha1*gradPhi/r(i)
      endif
      endif

C  if( (i==IA .or. i==IA/2) .and. (j==1 .or. j==JMax) ) then; -- for debug
C        write(*,'(a,3i5,1p4g12.3)')' irho0 i j C phi Fpsi H:',irho0,
C                            i,j,C, phi(i,j), psiOut,C -phi(i,j)-psiOut;
C  write(*,'(a,1p,4e20.12)')' on write delta, r^delta, theta :',
C                    delta, powrsmart(r(i), delta), theta(j),
C                     C -phi(i,j)-psiOut;
C       call mypause;
C  endif;
C*
Coutcylinder : write to channel 26+IB for plots and debug
C                              in cylinder coordinates  *
      enddo
      enddo
C-            Phirthet(0,0) = C -phi(0,0);
      write(*,*)'Phirthet 00 11:',Phirthet(0,0),Phirthet(1,1)
C-            call mypause;
C-          open(1,file='output',access='append')
C-          write(1,*) 'time',etime(tarray),dtime(tarray)
C- === PC ==>
C-          Call timer(ITime)
C-           IDTime=ITime-ITime0;
      t1=secnds(0.0)-t0
C- <== PC ===
C-           write(1,*) '*** IB, time', IB, IDTime;
C-            write(1,'(a,i6,1p,e11.3)') '#*** IB, time:', IB, t1;
      Re=sqrt(phiu/(G_N*rho0))
      if(Lrho)then
      rhokeep(irho0)=rho0
      Jkeep(irho0)=CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50
      Mkeep(irho0)=CM*Re**3*rho0/S_mass
      endif
      write(*,'(a,1p,e12.5)') 'Re=',Re;
      write(*,'(a,1p,e12.5)') 'phiu=',phiu;
      write(*,'(a,1p,e12.5)') 'G_N=',G_N;
      write(*,'(a,1p,e12.5)') 'rho0=',rho0;
      call mypause;
C --- WRITING IN ARRAY TTT (FOR MAIN PROGRAM). 14.06.2018.
      TTT(1)=rho01T         ! rho0
	TTT(2)=r(IB)**delta   ! oblatness
	TTT(3)=Re
	TTT(4)=CM*Re**3*rho0/S_mass  ! M_S
     	TTT(5)=CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50 ! J_CGS
      TTT(6)=T/abs(W)       ! TAU
	TTT(7)=Pip/Eth+one    ! gammaT
	TTT(8)=rhomax
	TTT(9)=phi(I0,JMax)
      TTT(10)=phi(IA,JMax)
	TTT(11)=phi(IB,JB)
	TTT(12)=CM
	TTT(13)=V
	TTT(14)=CJ
	TTT(15)=CI
	TTT(16)=W
	TTT(17)=VT
	TTT(18)=Cpsi

      print*, INDEXT
!	READ*
       IF (IOBHOD.EQ.1) GO TO 777
        write(1,22)INDEXT,rho01T,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	CJ*sqrt(G_N)
     *	*Re**5*rho0**1.5d0/1d50,T/abs(W),Pip/Eth+one,rhomax,phi(I0,JM
     *ax),phi(IA,JMax),phi(IB,JB),CM,V,CJ,CI,W,VT,Cpsi,IB
  777 CONTINUE

C-,t1;
C-          close(1)
      write(*,*)' output integrals'
      write(2,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     * ,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
C-            close(2);
      write(*,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      write(*,*)'CJ, VT',CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,VT
C-            pause;
C: prepare spline for changing to cylindrical coords *
C------- '-->Entering Node %_OutputResults_Prepspline:'
C-        NrDAT=NrDAT;
C-        NzDAT=NzDAT;
      write(*,*)' N=',NrDAT,'   NzDAT=',NzDAT,'   NtDAT=',NtDAT
      write(*,*)' NrDAT=',NrDAT,'   NzDAT=',NzDAT
C-        call mypause;
      do ir=0,NrDAT
      rspl(ir)=powrsmart(r(ir),delta)
      enddo
      do jt=0,NtDAT
      thet(jt)=min(theta(jt),pi/2.d0)
      enddo
C-     .. Executable Statements ..
      WRITE(*,*)'E01DAF Program Results'
C-     the number of X points, MX, and the values of the
C-     X co-ordinates.
C-        MX=NrDAT;
      MX=NrDAT+1
C-     the number of Y points, MY, and the values of the
C-     Y co-ordinates.
      MY=NtDAT+1
C-     Read the function values at the grid points.
C       do jt=0,MY;
C          do ir=0,MX;
C            Fspl((MY+1)*ir+jt+1)=Phirthet(ir,jt);
C          enddo;
C       enddo;*
      do jt=1,MY
      do ir=1,MX
      Fspl(MY*(ir-1)+jt)=Phirthet(ir-1,jt-1)
      enddo
      enddo
      IFAIL=0
C-
C-  *     Generate the (X,Y,F) interpolating bicubic B-spline.
C-       CALL E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
C-     Generate the (rspl,thet,Fspl) interpolating bicubic B-spline.
      CALL E01DAF(MX,MY,rspl,thet,Fspl,PX,PY,LAMDA,MU,Cspl,WRK,IFAIL)
C- i.e. here the function Fspl(rspl,thet) is given on a grid NrDAT+1 \times NtDAT+1
C- as an array, the knots are in LAMDA and MU 1D arrays:
C
C      SUBROUTINE E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC     Derived from DASL routine B2IRE.
CC     .. Parameters ..
C      CHARACTER*6       SRNAME
C      PARAMETER         (SRNAME='E01DAF')
C      DOUBLE PRECISION  ONE
C      PARAMETER         (ONE=1.0D+0)
CC     .. Scalar Arguments ..
C      INTEGER        IFAIL, MX, MY,
C         PX, PY -- output
CF(MX*MY) – real array  -- Input
COn entry: F(my * (q - 1)+r) must contain f_{q,r}, for q = 1, 2, ... ; mx , r = 1, 2, ... , my.
C
CC(MX*MY) – real array -- Output
COn exit: the coefﬁcients of the spline interpolant. C(my * (i - 1) + j)
Ccontains the coefﬁcient c_{ij} described in Section 3.
C
C
CC     .. Array Arguments ..
C      DOUBLE PRECISION  C(MX*MY), F(MX*MY), LAMDA(MX+4), MU(MY+4),
C     *                  WRK((MX+6)*(MY+6)), X(MX), Y(MY)
C_________________________________________________________________________
C*     E01DAF Example Program Text
C*     Mark 14 Release.  NAG Copyright 1989.
C*     .. Parameters ..
C      INTEGER          NIN, NOUT
C      PARAMETER        (NIN=5,NOUT=6)
C      INTEGER          MXMAX, MYMAX
C      PARAMETER        (MXMAX=20,MYMAX=MXMAX)
C      INTEGER          LIWRK, LWRK
C      PARAMETER        (LIWRK=MXMAX+2*(MXMAX-3)*(MYMAX-3),LWRK=(MXMAX+6)
C     +                 *(MYMAX+6))
C*     .. Local Scalars ..
C      DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO
C      INTEGER          I, IFAIL, J, MX, MY, NX, NY, PX, PY
C*     .. Local Arrays ..
C      DOUBLE PRECISION C(MXMAX*MYMAX), F(MXMAX*MYMAX), FG(MXMAX*MYMAX),
C     +                 LAMDA(MXMAX+4), MU(MYMAX+4), TX(MXMAX),
C     +                 TY(MYMAX), WRK(LWRK), X(MXMAX), Y(MYMAX)
C      INTEGER          IWRK(LIWRK)
C      CHARACTER*10     CLABS(MYMAX), RLABS(MXMAX)
C  *
C-
C --     Print the knot sets, LAMDA and MU.
C       WRITE (*,*)'               I   Knot LAMDA(I)      J     Knot MU(J)';
C       do  J = 4, MAX(PX,PY) - 3;
C          IF (J.LE.PX-3 .AND. J.LE.PY-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J), J, MU(J);
C          ELSE IF (J.LE.PX-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J);
C          ELSE IF (J.LE.PY-3) THEN;
C             WRITE (*, 99996) J, MU(J);
C          END IF;
C       enddo;*
C-        call mypause;
C-     Print the spline coefficients.
C-        WRITE (*,*) 'The B-Spline coefficients computed in Cspl';
C-        WRITE (*,99999) (Cspl(I),I=1,MX*MY);
C------- '<--Leaving  Node %_OutputResults_Prepspline:'
C:  FROM Phirthet interpolate to phiRcylZ on a grid of Rcylindrical and Z *
C------- '-->Entering Node %_OutputResults_RcylZ:'
C- X(1)=0.; -- is Rcyl coord
C now in %_init:
C   do icyl=1,NrDAT;
C     Rcyld(icyl)=rspl(icyl);
C--      X(i)=Rcyld(i); -- X(i) not needed
C   enddo;
C   Rcyld(0)=0.d0;
C*
C- Y(1)=0.; -- is Z coord
      Zd(0)=0.d0
      do jcyl=1,NzDAT
      Zd(jcyl)=(Rcyld(NrDAT)/real(NzDAT))*real(jcyl)
      write(*,*)' jcyl,  Zd(jcyl)=',jcyl,Zd(jcyl)
C-       Y(jcyl)=Zd(jcyl); -- Y(i) not needed
      enddo
      do icyl=0,NrDAT
      do jcyl=0,NzDAT
C- height loop
      radius=sqrt(Rcyld(icyl)**2+Zd(jcyl)**2)
      if(jcyl.GT.0.d0)then
      thetspl=atan(Rcyld(icyl)/Zd(jcyl))
      else
      thetspl=pi/2.d0
      endif
C-       write(*,'(a,2i4,a,1p,4e11.3)')' icyl jcyl=',icyl,jcyl,' X Y radius thetspl =',X(icyl),
C-          Y(jcyl),radius,thetspl;
Clin : find Phi for point icyl,jcyl linearly interpolating from
C                 Phirthet at radius,thetspl *
      Xrad(jcyl)=max(rspl(0),min(radius,rspl(NrDAT)))
      Ythet(jcyl)=max(thet(0),min(thetspl,thet(NtDAT)))
C- call mypause;
      enddo
C- jcyl, i.e. height loop
C-      Xrad(0)=0.d0;
C-      Ythet(0)=0.d0;
C: find PhiRcylZ(icyl,jcyl) for NrDAT points icyl, and current jcyl by spline interpolation from
C                 Phirthet at radius,thetspl *
C-
C-        Evaluate the spline at M points
      M=NzDAT+1
      CALL E02DEF(M,PX,PY,Xrad(0),Ythet(0),LAMDA,MU,Cspl,FF,WRK,IWRK,
     *	IFAIL)
C
C      SUBROUTINE E02DEF(M,PX,PY,X,Y,LAMDA,MU,Cspl,FF,WRK,IWRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC
CC     Derived from DASL routine B2VRE.
CC
CC     E02DEF. An algorithm for evaluating a bicubic polynomial
CC     spline S(X,Y) from its B-spline representation at the M
CC     points (X(I),Y(I)), I = 1, 2, ..., M.
CC
CC     Input Parameters:
CC        M        The number of evaluation points.
CC        PX       NXKNTS + 8,  where  NXKNTS  is number
CC                    of interior  X-knots.
CC        PY       NYKNTS + 8,  where  NYKNTS  is number
CC                    of interior  Y-knots.
CC        X        X-values.
CC        Y        Y-values.
CC        LAMDA    The X-knots.
CC        MU       The Y-knots.
CC        Cspl        B-spline coefficients of  S.  First
CC                    subscript relates to  X.
CC
CC     Output parameter:
CC        FF       Values of spline.
CC                 On exit, FF(I) contains the value of
CC                 the spline evaluated at point (X(I),Y(I)),
CC                 for I = 1,..,M.
CC
CC     Workspace parameters:
CC        WRK      Real workspace of dimension at least (PY-4).
CC        IWRK     Integer workspace of dimension at least (PY-4).
CC
CC     Failure indicator parameter:
CC        IFAIL    Failure indicator:
CC                 1 -  PX     .LT. 8,    or
CC                      PY     .LT. 8,    or
CC                      M      .LT. 1.
CC                 2 -  E02DFW  failure for  X  or  Y.
CC                 3 -  At least one point (X(K),Y(K)) lies
CC                      outside the rectangle defined by
CC                      LAMDA(4), LAMDA(PX-3), MU(4) and MU(PY-3).
C  *
      If(Ifail.ne.0)then
      write(*,*)' Ifail after spline E02DEF=',Ifail
      call mypause
      endif
      do jcyl=1,NzDAT+1
      phiRcylZ(icyl,jcyl-1)=FF(jcyl)
C- this is dp either phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      if(icyl.GT.1.and.jcyl.GT.1)phiRcylZplot(icyl-1,jcyl-1)=FF(jcyl)
C- this is real*4 phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      enddo
      enddo
C- icyl, i.e. rcyl loop
C------- '<--Leaving  Node %_OutputResults_RcylZ:'
C: plots and/or control integrals in cylindrical coords *
C------- '-->Entering Node %_OutputResults_controls:'
C -- for supermongo
C        open(3,file='PhiRcylZ.res',form='unformatted'); -- this file is for sm image
C--  now write your data into arr
C        write(3) NrDAT,NzDAT;
C        write(3) phiRcylZ;
C*
C- -- for pgplot - pgxtal
C-      CALL PLOT(X,NrDAT,Y,NzDAT,PhiRcylZ,NrDAT,NzDAT,Wdummy,1.5,4,'Rcyl-axis','Z-axis',
C-                'potential phi(Rcyl,Z)');
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
      dZcyl=Zd(2)-Zd(1)
C- for uniform Zcyl grid
      write(*,*)' Zd(0)=',Zd(0)
      write(*,*)' dZcyl,Zd(1)-Zd(0)=',dZcyl,Zd(1)-Zd(0)
C-        call mypause;
C: for Zcyl=Zd(0) compare rho and phi in spherical
C                 and cylinder coordinates in the same points*
C------- '-->Entering Node %_OutputResults_controls_equator:'
      write(*,*)' NrDAT,IA=',NrDAT,IA
C: try different 1D interpolation procedures from NAG to produce phiEqu(Rcyl)
C      for z=0 and plot them online with supermongo  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate:'
      write(*,*)' before OutputResults_controls_equator_interpolate
     *	_Chebyshev'
C-       call mypause;
C: use E01AEF to build Chebyshev interpolant and E02AKF to find interpolated values
C                  with  Chebyshev polynomials
C                  1) along z=0 for intermediate points of Rcyl grid;
C                  2) use phi on r,theta grid for constant r and interpolate phi(theta)  to theta=pi/2  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C-      E01AEF. A ROUTINE, WITH CHECKS, WHICH DETERMINES AND
C-      REFINES A POLYNOMIAL INTERPOLANT  Q(Xch)  TO DATA WHICH
C-      MAY CONTAIN DERIVATIVES.
C-
C-      INPUT PARAMETERS
C-         Mcheb        NUMBER OF DISTINCT Xch-VALUES
C-         XMIN,
C-         XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C-         Xch        INDEPENDENT VARIABLE VALUES (DISTINCT)
C-         Ych        VALUES AND DERIVATIVES OF
C-                     DEPENDENT VARIABLE.
C-         nDers       HIGHEST ORDER OF DERIVATIVE AT EACH Xch-VALUE.
C-         in the Rotat code we have no derivatives, so all nDers==0
C-         Ncheb        NUMBER OF INTERPOLATING CONDITIONS.
C-                     Ncheb == Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb).
C-         in the Rotat code we have no derivatives, so Ncheb == Mcheb
C-         ITMIN,
C-         ITMAX    MINIMUM AND MAXIMUM NUMBER OF ITERATIONS TO BE
C-                     PERFORMED.
C-
C-      OUTPUT PARAMETERS
C-         chebA        CHEBYSHEV COEFFICIENTS OF  Q(Xch)
C-
C-      WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C-         ChWRK      REAL WORKSPACE ARRAY.  THE FIRST IMAX ELEMENTS
C-                     CONTAIN, ON EXIT, PERFORMANCE INDICES FOR
C-                     THE INTERPOLATING POLYNOMIAL, AND THE NEXT
C-                     Ncheb  ELEMENTS THE COMPUTED RESIDUALS
C-         ChLWRK     DIMENSION OF ChWRK. ChLWRK MUST BE AT LEAST
C-                     7*Ncheb + 5*IMAX + Mcheb + 2, WHERE
C-                     IMAX IS ONE MORE THAN THE LARGEST ELEMENT
C-                     OF THE ARRAY nDers.
C-         IWRK     INTEGER WORKSPACE ARRAY.  ON EXIT,  IWRK(1)
C-                     CONTAINS THE NUMBER OF ITERATIONS TAKEN
C-         ChLIWRK    DIMENSION OF IWRK.  AT LEAST 2*Mcheb + 2.
C-
C-      FAILURE INDICATOR PARAMETER
C-         IFAIL    FAILURE INDICATOR.
C-                     0 - SUCCESSFUL TERMINATION.
C-                     1 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            Mcheb AT LEAST 1,
C-                            Ncheb = Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb),
C-                            ChLWRK AT LEAST 7*Ncheb + 5*IMAX + Mcheb + 2,
C-                            ChLIWRK AT LEAST 2*Mcheb + 2.
C-                     2 - FOR SOME I, nDers(I) IS LESS THAN 0.
C-                     3 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            XMIN STRICTLY LESS THAN XMAX,
C-                            FOR EACH I, Xch(I) MUST LIE IN THE
C-                               INTERVAL XMIN TO XMAX,
C-                            THE Xch-VALUES MUST ALL BE DISTINCT
C-                     4 - NOT ALL PERFORMANCE INDICES LESS THAN
C-                         ONE, BUT ITMAX ITERATIONS PERFORMED,
C-                     5 - COMPUTATION TERMINATED BECAUSE
C-                         ITERATIONS DIVERGING.
C-
C-
C-      CHECK AND SET ITERATION LIMITS
C-
C-      .. Parameters ..
C-      CHARACTER*6       SRNAME
C-      PARAMETER         (SRNAME='E01AEF')
C-      .. Scalar Arguments ..
C-      DOUBLE PRECISION  XMAX, XMIN
C-      INTEGER           IFAIL, ITMAX, ITMIN, ChLIWRK, ChLWRK, Mcheb, Ncheb
C-      .. Array Arguments ..
C-      DOUBLE PRECISION  chebA(Ncheb), WRK(ChLWRK), Xch(Mcheb), Ych(Ncheb)
C-      INTEGER           nDers(Mcheb), IWRK(ChLIWRK)
      write(*,*)' entering OutputResults_controls_equator_
     *	interpolate_Chebyshev'
C: define for each I from 1 to Mcheb the values of nDers(I), Xch(I),
C                and (Ych(J),J=N+1,N+nDers(I)+1).
C            nDersMAX = MAX(nDersMAX,nDers(I))
C            N = N + nDers(I) + 1
C   *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
      enddo
      Mcheb=NrDAT+1
C- too many points
      Mcheb=8
C- try to do in pieces (may be overlapping), a loop over  pieces will be needed
      skipCheb=20
      Mcheb=(NrDAT+1)/skipCheb
      write(*,*)'  NrDAT+1,skipCheb,Mcheb=',NrDAT+1,skipCheb,Mcheb
C-          pause;
      Ncheb=0
      do I=1,Mcheb
      nDers(I)=0
C- non-zero if derivatives are given
      Xch(I)=Rcyld(skipCheb*I-1)
      do J=Ncheb+1,Ncheb+nDers(I)+1
      Ych(J)=phiRcylZ(skipCheb*J-1,0)
      enddo
      Ncheb=Ncheb+nDers(I)+1
      enddo
      XMIN=Xch(1)
      XMAX=Xch(Mcheb)
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDerMX)THEN
      IFAIL=1
      nDersMAX=MAXVAL(nDers)
C- here must be zero
      write(*,*)' nDersMAX=',nDersMAX
C-             call mypause;
      ENDIF
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      write(*,*)' before   IF (Ncheb.LE.ChNMAX  ...)'
      write(*,*)' Ncheb.LE.ChNMAX .AND. nDersMAX.LE.nDersMX',Ncheb,
     *	ChNMAX,nDersMAX,nDersMX
C-          pause;
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDersMX)THEN
      IFAIL=1
      ITMIN=0
      ITMAX=0
C- The next call to E01AEF can be used for obtaining  interpolant q(x) as a series of Chebyshev
C- polynomials with coefficients chebA(i).
C- The polynomial interpolant can subsequently be evaluated for any value of x in the
C- given range by using
C- E02AKF. Chebyshev-series representations of the derivative(s) and integral(s) of q(x) may be
C- obtained by (repeated) use of E02AHF and E02AJF.
      CALL E01AEF(Mcheb,XMIN,XMAX,Xch,Ych,nDers,Ncheb,ITMIN,ITMAX,chebA
     *	,ChWRK,ChLWRK,IWRK,ChLIWRK,IFAIL)
      write(*,*)' after call E01AEF(Mcheb, ...)'
C-             pause;
      WRITE(NOUT,*)
      IF(IFAIL.EQ.0.OR.IFAIL.GE.4)THEN
      WRITE(NOUT,89999)'Total number of interpolating conditions ='
     *	,Ncheb
      WRITE(NOUT,*)
      WRITE(NOUT,*)'Interpolating polynomial'
      WRITE(NOUT,*)
      WRITE(NOUT,*)'   I    Chebyshev Coefficient chebA(I+1)'
      DO I=1,Ncheb
      WRITE(NOUT,89998)I-1,chebA(I)
      ENDDO
      WRITE(NOUT,*)
      WRITE(NOUT,*)'  Xch    R   Rth derivative    Residual'
      IY=0
      IRES=nDersMAX+1
      DO I=1,Mcheb
      nDers1=nDers(I)+1
      DO J=1,nDers1
      IY=IY+1
      IRES=IRES+1
      IF(J-1.NE.0)THEN
      WRITE(NOUT,89997)J-1,Ych(IY),ChWRK(IRES)
      ELSE
      WRITE(NOUT,89996)Xch(I),'   0',Ych(IY),ChWRK(IRES)
      ENDIF
      ENDDO
      ENDDO
      ELSE
      WRITE(NOUT,89995)'E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,Ncheb
     *	,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
      WRITE(NOUT,89995)'Chebyshev E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
      do i=1,MCheb
      IFAIL=0
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
C- this computes a function froom Chebyshev expansion:
      CALL E02AKF(NCheb+1,XMIN,XMAX,ChebA,1,NCheb+1,Xch(I),
     *	 resCheb(I),IFAIL)
      WRITE(*,'(a,i5,1p,2g12.3)')' ifail Xch(i) testCheb =',ifail,
     *	Xch(i),resCheb(i)
C- corresponds to ptype 6 3 in interactive sm:
      pt(i)=63.d0
      enddo
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
C-       CALL E01BAF(N,Xch,Ych,LAMDA,C,LCK,WRK,LWRK,IFAIL);
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C: for similar tasks E01BAF, or try the ready 2D spline *
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate:'
      Zcyl=Zd(0)
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      rho=EOS(phiu*(C-phia-psi),1)
C- in phys. units from EOS
C-            write(*,'(a,1p,4g15.5)')' Rcyl, r, phia, psi=',Rcyl, r(icyl), phia, psi;
C-            write(*,'(a,1p,g15.5)')' rho on equ=',rho;
C-            write(*,'(a,1p,g15.5)')' rho saved =',rhoOn_r_thetaGrid(icyl,JMax);
C-       call mypause;
      enddo
C------- '<--Leaving  Node %_OutputResults_controls_equator:'
C: output CMcyl, CJcyl etc for rho centered *
C------- '-->Entering Node %_OutputResults_controls_rhocentered:'
      CMcyl(0)=0.d0
C- mass within Rcyl(icyl)
      CJcyl(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(7,'(3a)')'#     Rcyla,              CMcyl(icyl),
     *	sigma','             sigmaGM,              ERROR','       IFA
     *IL,    CJcyl(icyl)'
      do icyl=1,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      if(icyl.GT.1)then
C- now this "if  " is not needed
      Rcylm=Rcyld(icyl-1)
      else
      Rcylm=0.d0
      endif
      Rcyla=half*(Rcyl+Rcylm)
C-            psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)));
      psi=Fpsi(Rcyla,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyla,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      CMcyl(icyl)=CMcyl(icyl-1)
      CJcyl(icyl)=CJcyl(icyl-1)
      sigma=0.d0
      do jcyl=1,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=(phiRcylZ(icyl-1,jcyl-1)+phiRcylZ(icyl-1,jcyl)+
     *	phiRcylZ(icyl,jcyl-1)+phiRcylZ(icyl,jcyl))/four
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dV=2.d0*(Rcyl**2-Rcylm**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      	sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      CMcyl(icyl)=CMcyl(icyl)+rho*dV
      CJcyl(icyl)=CJcyl(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd(1),rhoRcylZ(1,icyl),NzDAT,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyla
      write(7,'(1p,5g20.8,i5,g20.8)')Rcyla,CMcyl(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcyl(icyl)
      enddo
C- loop on Rcyl
      Call GillMiller(Rcyld(1),sigma2pir(1),NrDAT,CMGM,ERROR,IFAIL)
      write(7,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,
     *	ERROR,IFAIL
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(2,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      if(Lpolytrope)then
      callpolytrope(CN)
      write(2,2)'poly',rho0,r(IB)**delta/alphaPolytrope,Re/
     *	alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/
     *    alphaPolytrope)*
     **3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI*Re**5/(alphaPolytrope**5
     *),T/abs(W),Pip/Eth+one,VT
      endif
C-        pause;
C------- '<--Leaving  Node %_OutputResults_controls_rhocentered:'
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
C: output CMcylg, CJcylg etc for rho on cylindrical grid  *
C------- '-->Entering Node %_OutputResults_controls_rhoOnCylGrid:'
      CMcylg(0)=0.d0
C- mass within Rcyl(icyl)
      CJcylg(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(8,'(3a)')'#     Rcyla,              CMcylg(icyl),
     *   sigma','             sigmaGM,              ERROR','       IF
     *AIL,    CJcylg(icyl)'
      Rcylm=0.d0
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyl,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      if(icyl.GT.0)then
      CMcylg(icyl)=CMcylg(icyl-1)
      CJcylg(icyl)=CJcylg(icyl-1)
      Rcylm=Rcyld(icyl-1)
      endif
      sigma=0.d0
      dojcyl=0,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=phiRcylZ(icyl,jcyl)
C- here phia on cylindrical grid
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dRcyl=Rcyl-Rcylm
      if(jcyl.LT.NzDAT)then
      dV=2.d0*((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      else
      dV=((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- not multiplied by 2 since dZcyl is on border
      sigma=sigma+rho*dZcyl
C- not multiplied by 2
      endif
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      CMcylg(icyl)=CMcylg(icyl)+rho*dV
      CJcylg(icyl)=CJcylg(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd,rhoRcylZ(0,icyl),NzDAT+1,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyl
      write(8,'(1p,5g20.8,i5,g20.8)')Rcyl,CMcylg(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcylg(icyl)
      enddo
C- loop on Rcyl
      Call GillMiller(Rcyld,sigma2pir,NrDAT+1,CMGM,ERROR,IFAIL)
      write(8,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,
     *	ERROR,IFAIL
      do icyl=0,NrDAT
      mcyl(icyl)=CMcylg(icyl)/CMcylg(NrDAT)
      enddo
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(8,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     * 	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     *	V*Re**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),
     *    Pip/Eth+one,VT
      if(Lpolytrope)write(8,2)'poly',rho0,r(IB)**delta/alphaPolytrope,
     *	Re/alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/al
     *phaPolytrope)**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI*Re**5
     */(alphaPolytrope**5),
C- *rho0),
     *T/abs(W),Pip/Eth+one,VT
C-        pause;
      LinitMcyl=.false.
C- needed for jmscf -- 1st model is rigid, 2nd jmscf
C------- '<--Leaving  Node %_OutputResults_controls_rhoOnCylGrid:'
C------- '<--Leaving  Node %_OutputResults_controls:'
C------- '<--Leaving  Node %_OutputResults:'
      close(26+INDEXT)
      enddo
C-            write(4,@ap(1p,@irho2max g15.5)@ap)((rhokeep(krho),rhokeep(krho)),krho=0,irho0max);
C-            write(4,@ap(1p,@irho2max g15.5)@ap)((Jkeep(krho),Mkeep(krho)),krho=0,irho0max);
      write(4,' (1p,102g15.5)' )(Jkeep(krho),Mkeep(krho),krho=0,
     *	irho0max)
Cmine      elseif(rotlaw.EQ.'rigid')then
      elseif(rotlaw.eq.'rigid'.or.rotlaw.eq.'collap'.or.
     *	rotlaw.eq.'omega'.or.rotlaw.eq.'saved')then
C: read IB from tableIB as in Hachisu or in our paper
C                       AA290(1994)674 *
C------- '-->Entering Node %_tableIB:'
      do itable=1,NtableIB
      IB=nint(IA*tableIB(itable))
	print*,'start new cycle, OK?'
	read*
C- assuming IB<=9999
      IF(.NOT.(IB.LT.10))GOTO 09991
      write(str,'(a,I1)')'000',IB
      GOTO 09988
09991 CONTINUE
      IF(.NOT.(IB.LT.100))GOTO 09990
      write(str,'(a,I2)')'00',IB
      GOTO 09988
09990 CONTINUE
      IF(.NOT.(IB.LT.1000))GOTO 09989
      write(str,'(a,I3)')'0',IB
      GOTO 09988
09989 CONTINUE
      write(str,'(I4)')IB
09988 CONTINUE
      open(26+IB,file='phi'//str(1:len_trim(str))//'.res')
C: iterate phi etc. for rigid rotation *
C------- '-->Entering Node %_FindPhi:'
      write(*,*)' entering _FindPhi'
C-            call mypause;
      LogInit=.true.
      write(*,'(a,1p,g15.5,i7)')'rho0, IB',rho0,IB
      write(*,'(a,1p,3g15.5)')'phi(IA,JMax),phi(I0,JMax) =',
     *	phi(IA,JMax),phi(I0,JMax)
C-           'phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)))=',phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)));
C- write of Fpsi is forbidden if there is write in it
      write(*,'(2a)')' rotlaw=',rotlaw(1:len_trim(rotlaw))
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)
     *	*EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      write(*,'(a,1p,3g15.5)')'phiu,EOS(rho0,3),C=',phiu,EOS(rho0,3),C
C-            call mypause;
      n_max=300
!      epsMain=1.d-12
!      epsMain=1.d-9
        epsMain=1.d-8

C- to data file
      LogIter=.true.
      nStep=1
      do while(LogIter.and.(nStep.LE.n_max))
      C0=C
      Call new_phi(
C-r,alpha,theta,dtheta, -- now in module cylinderGrid
     *C,LogInit,LogRapid,rotlaw(1:len_trim(rotlaw)))
C- last argument is the rotation law
C              if(C>C0)then;
C                C=C0+min(abs(C-C0),0.03d0*abs(C0));
C              else;
C                C=C0-min(abs(C-C0),0.03d0*abs(C0));
C              endif;*
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18

!       pRINT*, '========2'
!	READ*

      if(Dabs((C0-C)/C).LT.epsMain)LogIter=.false.
      write(*,*)' nStep,C,phi(I0,JMax):',nStep,C,phi(I0,JMax)
      nStep=nStep+1
      enddo
C------- '<--Leaving  Node %_FindPhi:'
C: to channel 26+irho0 after convergence, find integrals,
C                                remap onto cylindrical coordinates,
C                                prepare plots etc. *
C------- '-->Entering Node %_OutputResults:'
      Call Integrals(CM,V,CJ,CI,T,W,Eth,Pip,C,omegaA,
     *	rotlaw(1:len_trim(rotlaw)))
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' VT=',VT
      rhomax=zero
      do i=0,IA
C- loop on r
      do j=0,JMax
C- loop on theta
      if(j.EQ.JMax)then
      phiEqu(i)=phi(i,j)
C                   write(*,'(a,2i5,1p,2g15.6)')' i, j, theta(j),phiEqu(i)=',
C                            i,j,theta(j),phiEqu(i);
C                   if(mod(i,20)==0)call mypause;*
      endif
C- write(*,*)i,j,delta,C;
      rho=EOS(phiu*(C-phi(i,j)-Fpsi(powrsmart(r(i),delta)*sin(theta(j))
     *	,rotlaw(1:len_trim(rotlaw)))),1)
C- in phys. units from EOS
      rhoOn_r_thetaGrid(i,j)=rho
C- save for future
      rhomax=max(rhomax,rho/rho0)
      if(i.GT.0.and.j.EQ.Jmax.and.irho0.EQ.irho0max)then
      H=EOS(rho,3)
      P=EOS(H,4)
      psi=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
      vphi2=-Fdpsi_drcyl(powrsmart(r(i),delta)*sin(theta(j)),
     *rotlaw(1:len_trim(rotlaw)))*powrsmart(r(i),delta)*sin(theta(JMax))
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,vphi2))
      write(12,'(1p,7g12.4)')r(i),rho,P,H,phi(i,j),psi,vphi
      endif
C: write to channel 26+irho0 for plots and debug
C                              in spherical coordinates  *
C- physical radius is r(i)**delta, i.e. powrsmart(r(i), delta)
      psiOut=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))
C-                  Phirthet(i,j) = C -phi(i,j)-psiOut;
      Phirthet(i,j)=phi(i,j)
      write(26+irho0,'(1p,3e20.12)')powrsmart(r(i),delta),theta(j),
     *	C-phi(i,j)-psiOut
C SB added
      if(j.EQ.JMax)then
      if(i.EQ.0)then
      write(13,'(a,i5)')'# IB=',IB
      write(13,'(6(a,9x))')'#  i','  r','alpha1','alpha2','gradPhi',
     *	'Omega^2'
      else
C                     write(*,*) 'i,r(i),theta(j)=',i,r(i),theta(j);
C                     write(*,'(a,2i5,1p4g12.3)')'i j phi Fpsi Fpsi/phi:',
C                            i,j,phi(i,j), psiOut, psiOut/(phi(i,j)-phi(I0,JMax));
C                            -- Phi0 == Phi_0 == phi(I0,JMax)
C                     write(*,'(a,1p4g12.3)')'phiB, phiA, Phi0, (phiB-phiA)/(PhiA-Phi0) :',
C                            phi(IB,0), phi(IA,JMax), phi(I0,JMax),
C                            (phi(IB,0)-phi(IA,JMax))/(phi(IA,JMax)-phi(I0,JMax));
C                     pause;
C                            *
      alpha1=psiOut/(phi(i,j)-phi(I0,JMax))
      alpha2=(phi(IB,0)-phi(IA,JMax))/(phi(IA,JMax)-phi(I0,JMax))
      gradPhi=(phi(i+1,Jmax)-phi(i-1,Jmax))/(2d0*dr)
      write(13,'(i5,1p,4g15.5,g17.5)')i,r(i),alpha1,alpha2,gradPhi,
     *	-alpha1*gradPhi/r(i)
      endif
      endif

C
C  if( (i==IA .or. i==IA/2) .and. (j==1 .or. j==JMax) ) then; -- for debug
C        write(*,'(a,3i5,1p4g12.3)')' irho0 i j C phi Fpsi H:',irho0,
C                            i,j,C, phi(i,j), psiOut,C -phi(i,j)-psiOut;
C  write(*,'(a,1p,4e20.12)')' on write delta, r^delta, theta :',
C                    delta, powrsmart(r(i), delta), theta(j),
C                     C -phi(i,j)-psiOut;
C       call mypause;
C  endif;
C*
Coutcylinder : write to channel 26+IB for plots and debug
C                              in cylinder coordinates  *
      enddo
      enddo
C-            Phirthet(0,0) = C -phi(0,0);
      write(*,*)'Phirthet 00 11:',Phirthet(0,0),Phirthet(1,1)
C-            call mypause;
C-          open(1,file='output',access='append')
C-          write(1,*) 'time',etime(tarray),dtime(tarray)
C- === PC ==>
C-          Call timer(ITime)
C-           IDTime=ITime-ITime0;
      t1=secnds(0.0)-t0
C- <== PC ===
C-           write(1,*) '*** IB, time', IB, IDTime;
C-            write(1,'(a,i6,1p,e11.3)') '#*** IB, time:', IB, t1;
      Re=sqrt(phiu/(G_N*rho0))
      if(Lrho)then
      rhokeep(irho0)=rho0
      Jkeep(irho0)=CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50
      Mkeep(irho0)=CM*Re**3*rho0/S_mass
      endif
      write(*,*)'Re=',Re
C-            call mypause;
      write(1,22)rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,CJ*sqrt(G_N)
     *	*Re**5*rho0**1.5d0/1d50,T/abs(W),Pip/Eth+one,rhomax,phi(I0,JM
     *ax),phi(IA,JMax),phi(IB,JB),CM,V,CJ,CI,W,VT,Cpsi,IB
C-,t1;
C-          close(1)
      write(*,*)' output integrals'
      write(2,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
C-            close(2);
      write(*,2)'intg',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      write(*,*)'CJ, VT',CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,VT
C-            pause;
C: prepare spline for changing to cylindrical coords *
C------- '-->Entering Node %_OutputResults_Prepspline:'
C-        NrDAT=NrDAT;
C-        NzDAT=NzDAT;
      write(*,*)' N=',NrDAT,'   NzDAT=',NzDAT,'   NtDAT=',NtDAT
      write(*,*)' NrDAT=',NrDAT,'   NzDAT=',NzDAT
C-        call mypause;
      do ir=0,NrDAT
      rspl(ir)=powrsmart(r(ir),delta)
      enddo
      do jt=0,NtDAT
      thet(jt)=min(theta(jt),pi/2.d0)
      enddo
C-     .. Executable Statements ..
      WRITE(*,*)'E01DAF Program Results'
C-     the number of X points, MX, and the values of the
C-     X co-ordinates.
C-        MX=NrDAT;
      MX=NrDAT+1
C-     the number of Y points, MY, and the values of the
C-     Y co-ordinates.
      MY=NtDAT+1
C-     Read the function values at the grid points.
C       do jt=0,MY;
C          do ir=0,MX;
C            Fspl((MY+1)*ir+jt+1)=Phirthet(ir,jt);
C          enddo;
C       enddo;*
      do jt=1,MY
      do ir=1,MX
      Fspl(MY*(ir-1)+jt)=Phirthet(ir-1,jt-1)
      enddo
      enddo
      IFAIL=0
C-
C-  *     Generate the (X,Y,F) interpolating bicubic B-spline.
C-       CALL E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
C-     Generate the (rspl,thet,Fspl) interpolating bicubic B-spline.
      CALL E01DAF(MX,MY,rspl,thet,Fspl,PX,PY,LAMDA,MU,Cspl,WRK,IFAIL)
C- i.e. here the function Fspl(rspl,thet) is given on a grid NrDAT+1 \times NtDAT+1
C- as an array, the knots are in LAMDA and MU 1D arrays:
C
C      SUBROUTINE E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC     Derived from DASL routine B2IRE.
CC     .. Parameters ..
C      CHARACTER*6       SRNAME
C      PARAMETER         (SRNAME='E01DAF')
C      DOUBLE PRECISION  ONE
C      PARAMETER         (ONE=1.0D+0)
CC     .. Scalar Arguments ..
C      INTEGER        IFAIL, MX, MY,
C         PX, PY -- output
CF(MX*MY) – real array  -- Input
COn entry: F(my * (q - 1)+r) must contain f_{q,r}, for q = 1, 2, ... ; mx , r = 1, 2, ... , my.
C
CC(MX*MY) – real array -- Output
COn exit: the coefﬁcients of the spline interpolant. C(my * (i - 1) + j)
Ccontains the coefﬁcient c_{ij} described in Section 3.
C
C
CC     .. Array Arguments ..
C      DOUBLE PRECISION  C(MX*MY), F(MX*MY), LAMDA(MX+4), MU(MY+4),
C     *                  WRK((MX+6)*(MY+6)), X(MX), Y(MY)
C_________________________________________________________________________
C*     E01DAF Example Program Text
C*     Mark 14 Release.  NAG Copyright 1989.
C*     .. Parameters ..
C      INTEGER          NIN, NOUT
C      PARAMETER        (NIN=5,NOUT=6)
C      INTEGER          MXMAX, MYMAX
C      PARAMETER        (MXMAX=20,MYMAX=MXMAX)
C      INTEGER          LIWRK, LWRK
C      PARAMETER        (LIWRK=MXMAX+2*(MXMAX-3)*(MYMAX-3),LWRK=(MXMAX+6)
C     +                 *(MYMAX+6))
C*     .. Local Scalars ..
C      DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO
C      INTEGER          I, IFAIL, J, MX, MY, NX, NY, PX, PY
C*     .. Local Arrays ..
C      DOUBLE PRECISION C(MXMAX*MYMAX), F(MXMAX*MYMAX), FG(MXMAX*MYMAX),
C     +                 LAMDA(MXMAX+4), MU(MYMAX+4), TX(MXMAX),
C     +                 TY(MYMAX), WRK(LWRK), X(MXMAX), Y(MYMAX)
C      INTEGER          IWRK(LIWRK)
C      CHARACTER*10     CLABS(MYMAX), RLABS(MXMAX)
C  *
C-
C --     Print the knot sets, LAMDA and MU.
C       WRITE (*,*)'               I   Knot LAMDA(I)      J     Knot MU(J)';
C       do  J = 4, MAX(PX,PY) - 3;
C          IF (J.LE.PX-3 .AND. J.LE.PY-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J), J, MU(J);
C          ELSE IF (J.LE.PX-3) THEN;
C             WRITE (*, 99997) J, LAMDA(J);
C          ELSE IF (J.LE.PY-3) THEN;
C             WRITE (*, 99996) J, MU(J);
C          END IF;
C       enddo;*
C-        call mypause;
C-     Print the spline coefficients.
C-        WRITE (*,*) 'The B-Spline coefficients computed in Cspl';
C-        WRITE (*,99999) (Cspl(I),I=1,MX*MY);
C------- '<--Leaving  Node %_OutputResults_Prepspline:'
C:  FROM Phirthet interpolate to phiRcylZ on a grid of Rcylindrical and Z *
C------- '-->Entering Node %_OutputResults_RcylZ:'
C- X(1)=0.; -- is Rcyl coord
C now in %_init:
C   do icyl=1,NrDAT;
C     Rcyld(icyl)=rspl(icyl);
C--      X(i)=Rcyld(i); -- X(i) not needed
C   enddo;
C   Rcyld(0)=0.d0;
C*
C- Y(1)=0.; -- is Z coord
      Zd(0)=0.d0
      do jcyl=1,NzDAT
      Zd(jcyl)=(Rcyld(NrDAT)/real(NzDAT))*real(jcyl)
      write(*,*)' jcyl,  Zd(jcyl)=',jcyl,Zd(jcyl)
C-       Y(jcyl)=Zd(jcyl); -- Y(i) not needed
      enddo
      do icyl=0,NrDAT
      do jcyl=0,NzDAT
C- height loop
      radius=sqrt(Rcyld(icyl)**2+Zd(jcyl)**2)
      if(jcyl.GT.0.d0)then
      thetspl=atan(Rcyld(icyl)/Zd(jcyl))
      else
      thetspl=pi/2.d0
      endif
C-       write(*,'(a,2i4,a,1p,4e11.3)')' icyl jcyl=',icyl,jcyl,' X Y radius thetspl =',X(icyl),
C-          Y(jcyl),radius,thetspl;
Clin : find Phi for point icyl,jcyl linearly interpolating from
C                 Phirthet at radius,thetspl *
      Xrad(jcyl)=max(rspl(0),min(radius,rspl(NrDAT)))
      Ythet(jcyl)=max(thet(0),min(thetspl,thet(NtDAT)))
C- call mypause;
      enddo
C- jcyl, i.e. height loop
C-      Xrad(0)=0.d0;
C-      Ythet(0)=0.d0;
C: find PhiRcylZ(icyl,jcyl) for NrDAT points icyl, and current jcyl by spline interpolation from
C                 Phirthet at radius,thetspl *
C-
C-        Evaluate the spline at M points
      M=NzDAT+1
      CALL E02DEF(M,PX,PY,Xrad(0),Ythet(0),LAMDA,MU,Cspl,FF,WRK,IWRK,
     *	IFAIL)
C
C      SUBROUTINE E02DEF(M,PX,PY,X,Y,LAMDA,MU,Cspl,FF,WRK,IWRK,IFAIL)
CC     MARK 14 RELEASE. NAG COPYRIGHT 1989.
CC
CC     Derived from DASL routine B2VRE.
CC
CC     E02DEF. An algorithm for evaluating a bicubic polynomial
CC     spline S(X,Y) from its B-spline representation at the M
CC     points (X(I),Y(I)), I = 1, 2, ..., M.
CC
CC     Input Parameters:
CC        M        The number of evaluation points.
CC        PX       NXKNTS + 8,  where  NXKNTS  is number
CC                    of interior  X-knots.
CC        PY       NYKNTS + 8,  where  NYKNTS  is number
CC                    of interior  Y-knots.
CC        X        X-values.
CC        Y        Y-values.
CC        LAMDA    The X-knots.
CC        MU       The Y-knots.
CC        Cspl        B-spline coefficients of  S.  First
CC                    subscript relates to  X.
CC
CC     Output parameter:
CC        FF       Values of spline.
CC                 On exit, FF(I) contains the value of
CC                 the spline evaluated at point (X(I),Y(I)),
CC                 for I = 1,..,M.
CC
CC     Workspace parameters:
CC        WRK      Real workspace of dimension at least (PY-4).
CC        IWRK     Integer workspace of dimension at least (PY-4).
CC
CC     Failure indicator parameter:
CC        IFAIL    Failure indicator:
CC                 1 -  PX     .LT. 8,    or
CC                      PY     .LT. 8,    or
CC                      M      .LT. 1.
CC                 2 -  E02DFW  failure for  X  or  Y.
CC                 3 -  At least one point (X(K),Y(K)) lies
CC                      outside the rectangle defined by
CC                      LAMDA(4), LAMDA(PX-3), MU(4) and MU(PY-3).
C  *
      If(Ifail.ne.0)then
      write(*,*)' Ifail after spline E02DEF=',Ifail
      call mypause
      endif
      do jcyl=1,NzDAT+1
      phiRcylZ(icyl,jcyl-1)=FF(jcyl)
C- this is dp either phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      if(icyl.GT.1.and.jcyl.GT.1)phiRcylZplot(icyl-1,jcyl-1)=FF(jcyl)
C- this is real*4 phi or C -phi(icyl,jcyl)-psiOut on Rcyl, Z grid
      enddo
      enddo
C- icyl, i.e. rcyl loop
C------- '<--Leaving  Node %_OutputResults_RcylZ:'
C: plots and/or control integrals in cylindrical coords *
C------- '-->Entering Node %_OutputResults_controls:'
C -- for supermongo
C        open(3,file='PhiRcylZ.res',form='unformatted'); -- this file is for sm image
C--  now write your data into arr
C        write(3) NrDAT,NzDAT;
C        write(3) phiRcylZ;
C*
C- -- for pgplot - pgxtal
C-      CALL PLOT(X,NrDAT,Y,NzDAT,PhiRcylZ,NrDAT,NzDAT,Wdummy,1.5,4,'Rcyl-axis','Z-axis',
C-                'potential phi(Rcyl,Z)');
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
      dZcyl=Zd(2)-Zd(1)
C- for uniform Zcyl grid
      write(*,*)' Zd(0)=',Zd(0)
      write(*,*)' dZcyl,Zd(1)-Zd(0)=',dZcyl,Zd(1)-Zd(0)
C-        call mypause;
C: for Zcyl=Zd(0) compare rho and phi in spherical
C                 and cylinder coordinates in the same points*
C------- '-->Entering Node %_OutputResults_controls_equator:'
      write(*,*)' NrDAT,IA=',NrDAT,IA
C: try different 1D interpolation procedures from NAG to produce phiEqu(Rcyl)
C      for z=0 and plot them online with supermongo  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate:'
      write(*,*)' before OutputResults_controls_equator_interpolate
     *	_Chebyshev'
C-       call mypause;
C: use E01AEF to build Chebyshev interpolant and E02AKF to find interpolated values
C                  with  Chebyshev polynomials
C                  1) along z=0 for intermediate points of Rcyl grid;
C                  2) use phi on r,theta grid for constant r and interpolate phi(theta)  to theta=pi/2  *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C-      E01AEF. A ROUTINE, WITH CHECKS, WHICH DETERMINES AND
C-      REFINES A POLYNOMIAL INTERPOLANT  Q(Xch)  TO DATA WHICH
C-      MAY CONTAIN DERIVATIVES.
C-
C-      INPUT PARAMETERS
C-         Mcheb        NUMBER OF DISTINCT Xch-VALUES
C-         XMIN,
C-         XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C-         Xch        INDEPENDENT VARIABLE VALUES (DISTINCT)
C-         Ych        VALUES AND DERIVATIVES OF
C-                     DEPENDENT VARIABLE.
C-         nDers       HIGHEST ORDER OF DERIVATIVE AT EACH Xch-VALUE.
C-         in the Rotat code we have no derivatives, so all nDers==0
C-         Ncheb        NUMBER OF INTERPOLATING CONDITIONS.
C-                     Ncheb == Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb).
C-         in the Rotat code we have no derivatives, so Ncheb == Mcheb
C-         ITMIN,
C-         ITMAX    MINIMUM AND MAXIMUM NUMBER OF ITERATIONS TO BE
C-                     PERFORMED.
C-
C-      OUTPUT PARAMETERS
C-         chebA        CHEBYSHEV COEFFICIENTS OF  Q(Xch)
C-
C-      WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C-         ChWRK      REAL WORKSPACE ARRAY.  THE FIRST IMAX ELEMENTS
C-                     CONTAIN, ON EXIT, PERFORMANCE INDICES FOR
C-                     THE INTERPOLATING POLYNOMIAL, AND THE NEXT
C-                     Ncheb  ELEMENTS THE COMPUTED RESIDUALS
C-         ChLWRK     DIMENSION OF ChWRK. ChLWRK MUST BE AT LEAST
C-                     7*Ncheb + 5*IMAX + Mcheb + 2, WHERE
C-                     IMAX IS ONE MORE THAN THE LARGEST ELEMENT
C-                     OF THE ARRAY nDers.
C-         IWRK     INTEGER WORKSPACE ARRAY.  ON EXIT,  IWRK(1)
C-                     CONTAINS THE NUMBER OF ITERATIONS TAKEN
C-         ChLIWRK    DIMENSION OF IWRK.  AT LEAST 2*Mcheb + 2.
C-
C-      FAILURE INDICATOR PARAMETER
C-         IFAIL    FAILURE INDICATOR.
C-                     0 - SUCCESSFUL TERMINATION.
C-                     1 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            Mcheb AT LEAST 1,
C-                            Ncheb = Mcheb + nDers(1) + nDers(2) + ... + nDers(Mcheb),
C-                            ChLWRK AT LEAST 7*Ncheb + 5*IMAX + Mcheb + 2,
C-                            ChLIWRK AT LEAST 2*Mcheb + 2.
C-                     2 - FOR SOME I, nDers(I) IS LESS THAN 0.
C-                     3 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C-                         HAS BEEN VIOLATED -
C-                            XMIN STRICTLY LESS THAN XMAX,
C-                            FOR EACH I, Xch(I) MUST LIE IN THE
C-                               INTERVAL XMIN TO XMAX,
C-                            THE Xch-VALUES MUST ALL BE DISTINCT
C-                     4 - NOT ALL PERFORMANCE INDICES LESS THAN
C-                         ONE, BUT ITMAX ITERATIONS PERFORMED,
C-                     5 - COMPUTATION TERMINATED BECAUSE
C-                         ITERATIONS DIVERGING.
C-
C-
C-      CHECK AND SET ITERATION LIMITS
C-
C-      .. Parameters ..
C-      CHARACTER*6       SRNAME
C-      PARAMETER         (SRNAME='E01AEF')
C-      .. Scalar Arguments ..
C-      DOUBLE PRECISION  XMAX, XMIN
C-      INTEGER           IFAIL, ITMAX, ITMIN, ChLIWRK, ChLWRK, Mcheb, Ncheb
C-      .. Array Arguments ..
C-      DOUBLE PRECISION  chebA(Ncheb), WRK(ChLWRK), Xch(Mcheb), Ych(Ncheb)
C-      INTEGER           nDers(Mcheb), IWRK(ChLIWRK)
      write(*,*)' entering OutputResults_controls_equator_
     *	interpolate_Chebyshev'
C: define for each I from 1 to Mcheb the values of nDers(I), Xch(I),
C                and (Ych(J),J=N+1,N+nDers(I)+1).
C            nDersMAX = MAX(nDersMAX,nDers(I))
C            N = N + nDers(I) + 1
C   *
C------- '-->Entering Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
      enddo
      Mcheb=NrDAT+1
C- too many points
      Mcheb=8
C- try to do in pieces (may be overlapping), a loop over  pieces will be needed
      skipCheb=20
      Mcheb=(NrDAT+1)/skipCheb
      write(*,*)'  NrDAT+1,skipCheb,Mcheb=',NrDAT+1,skipCheb,Mcheb
C-          pause;
      Ncheb=0
      do I=1,Mcheb
      nDers(I)=0
C- non-zero if derivatives are given
      Xch(I)=Rcyld(skipCheb*I-1)
      do J=Ncheb+1,Ncheb+nDers(I)+1
      Ych(J)=phiRcylZ(skipCheb*J-1,0)
      enddo
      Ncheb=Ncheb+nDers(I)+1
      enddo
      XMIN=Xch(1)
      XMAX=Xch(Mcheb)
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDerMX)THEN
      IFAIL=1
      nDersMAX=MAXVAL(nDers)
C- here must be zero
      write(*,*)' nDersMAX=',nDersMAX
C-             call mypause;
      ENDIF
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev_parChebyshev:'
      write(*,*)' before   IF (Ncheb.LE.ChNMAX  ...)'
      write(*,*)' Ncheb.LE.ChNMAX .AND. nDersMAX.LE.nDersMX',Ncheb,
     *	ChNMAX,nDersMAX,nDersMX
C-          pause;
      IF(Ncheb.LE.ChNMAX.AND.nDersMAX.LE.nDersMX)THEN
      IFAIL=1
      ITMIN=0
      ITMAX=0
C- The next call to E01AEF can be used for obtaining  interpolant q(x) as a series of Chebyshev
C- polynomials with coefficients chebA(i).
C- The polynomial interpolant can subsequently be evaluated for any value of x in the
C- given range by using
C- E02AKF. Chebyshev-series representations of the derivative(s) and integral(s) of q(x) may be
C- obtained by (repeated) use of E02AHF and E02AJF.
      CALL E01AEF(Mcheb,XMIN,XMAX,Xch,Ych,nDers,Ncheb,ITMIN,ITMAX,chebA,
     *	ChWRK,ChLWRK,IWRK,ChLIWRK,IFAIL)
      write(*,*)' after call E01AEF(Mcheb, ...)'
C-             pause;
      WRITE(NOUT,*)
      IF(IFAIL.EQ.0.OR.IFAIL.GE.4)THEN
      WRITE(NOUT,89999)'Total number of interpolating conditions =',
     *	Ncheb
      WRITE(NOUT,*)
      WRITE(NOUT,*)'Interpolating polynomial'
      WRITE(NOUT,*)
      WRITE(NOUT,*)'   I    Chebyshev Coefficient chebA(I+1)'
      DO I=1,Ncheb
      WRITE(NOUT,89998)I-1,chebA(I)
      ENDDO
      WRITE(NOUT,*)
      WRITE(NOUT,*)'  Xch    R   Rth derivative    Residual'
      IY=0
      IRES=nDersMAX+1
      DO I=1,Mcheb
      nDers1=nDers(I)+1
      DO J=1,nDers1
      IY=IY+1
      IRES=IRES+1
      IF(J-1.NE.0)THEN
      WRITE(NOUT,89997)J-1,Ych(IY),ChWRK(IRES)
      ELSE
      WRITE(NOUT,89996)Xch(I),'   0',Ych(IY),ChWRK(IRES)
      ENDIF
      ENDDO
      ENDDO
      ELSE
      WRITE(NOUT,89995)'E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
      WRITE(NOUT,89995)'Chebyshev E01AEF exits with IFAIL =',IFAIL
      write(*,*)' Ifail after Chebyshev E01AEF=',Ifail
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
      write(*,'(a,1p,(5g12.3))')' Xch(I)=',(Xch(I),I=1,MCheb)
      do i=1,MCheb
      IFAIL=0
      write(*,'(a,2i5,1p,2g12.3)')' Mcheb Ncheb Xmin Xmax =',Mcheb,
     *	Ncheb,XMIN,XMAX
C- this computes a function froom Chebyshev expansion:
      CALL E02AKF(NCheb+1,XMIN,XMAX,ChebA,1,NCheb+1,Xch(I),resCheb(I),
     *	IFAIL)
      WRITE(*,'(a,i5,1p,2g12.3)')' ifail Xch(i) testCheb =',ifail,
     *	Xch(i),resCheb(i)
C- corresponds to ptype 6 3 in interactive sm:
      pt(i)=63.d0
      enddo
C-                call mypause;
C-                pause; -- better than mypause, because exits smoothly on non-enter
      ENDIF
C-       CALL E01BAF(N,Xch,Ych,LAMDA,C,LCK,WRK,LWRK,IFAIL);
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate_Chebyshev:'
C: for similar tasks E01BAF, or try the ready 2D spline *
C------- '<--Leaving  Node %_OutputResults_controls_equator_interpolate:'
      Zcyl=Zd(0)
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
      phia=phiRcylZ(icyl,0)
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      rho=EOS(phiu*(C-phia-psi),1)
C- in phys. units from EOS
C-            write(*,'(a,1p,4g15.5)')' Rcyl, r, phia, psi=',Rcyl, r(icyl), phia, psi;
C-            write(*,'(a,1p,g15.5)')' rho on equ=',rho;
C-            write(*,'(a,1p,g15.5)')' rho saved =',rhoOn_r_thetaGrid(icyl,JMax);
C-       call mypause;
      enddo
C------- '<--Leaving  Node %_OutputResults_controls_equator:'
C: output CMcyl, CJcyl etc for rho centered *
C------- '-->Entering Node %_OutputResults_controls_rhocentered:'
      CMcyl(0)=0.d0
C- mass within Rcyl(icyl)
      CJcyl(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(7,'(3a)')'#     Rcyla,              CMcyl(icyl),
     *	sigma','             sigmaGM,              ERROR','       IFA
     *IL,    CJcyl(icyl)'
      do icyl=1,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      if(icyl.GT.1)then
C- now this "if  " is not needed
      Rcylm=Rcyld(icyl-1)
      else
      Rcylm=0.d0
      endif
      Rcyla=half*(Rcyl+Rcylm)
C-            psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)));
      psi=Fpsi(Rcyla,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyla,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      CMcyl(icyl)=CMcyl(icyl-1)
      CJcyl(icyl)=CJcyl(icyl-1)
      sigma=0.d0
      do jcyl=1,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=(phiRcylZ(icyl-1,jcyl-1)+phiRcylZ(icyl-1,jcyl)+
     *	phiRcylZ(icyl,jcyl-1)+phiRcylZ(icyl,jcyl))/four
C-                rho=EOS(phiu*(phiRcylZ(icyl,jcyl)),1)/rho0;
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dV=2.d0*(Rcyl**2-Rcylm**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      	sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      CMcyl(icyl)=CMcyl(icyl)+rho*dV
      CJcyl(icyl)=CJcyl(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd(1),rhoRcylZ(1,icyl),NzDAT,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyla
      write(7,'(1p,5g20.8,i5,g20.8)')Rcyla,CMcyl(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcyl(icyl)
      enddo
C- loop on Rcyl
      CallGillMiller(Rcyld(1),sigma2pir(1),NrDAT,CMGM,ERROR,IFAIL)
      write(7,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,
     *	ERROR,IFAIL
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(2,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one,VT
      if(Lpolytrope)then
      call polytrope(CN)
      write(2,2)'poly',rho0,r(IB)**delta/alphaPolytrope,Re/
     *	alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/
     *    alphaPolytrope)**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,
     *    CI*Re**5/(alphaPolytrope**5),
C- *rho0),
     *T/abs(W),Pip/Eth+one,VT
      endif
C-        pause;
C------- '<--Leaving  Node %_OutputResults_controls_rhocentered:'
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
C: output CMcylg, CJcylg etc for rho on cylindrical grid  *
C------- '-->Entering Node %_OutputResults_controls_rhoOnCylGrid:'
      CMcylg(0)=0.d0
C- mass within Rcyl(icyl)
      CJcylg(0)=0.d0
C- angular momentum within Rcyl(icyl)
      write(8,'(3a)')'#     Rcyla,              CMcylg(icyl),
     *       sigma','             sigmaGM,              ERROR','
     *     IF AIL,    CJcylg(icyl)'
      Rcylm=0.d0
      do icyl=0,NrDAT
      Rcyl=Rcyld(icyl)
C- here dp
      psi=Fpsi(Rcyl,rotlaw(1:len_trim(rotlaw)))
      c1=-Fdpsi_drcyl(Rcyl,rotlaw(1:len_trim(rotlaw)))*Rcyla
C- it is minus, now the sign is opposite in chi -- check
      vphi=sqrt(max(zero,c1))
      if(icyl.GT.0)then
      CMcylg(icyl)=CMcylg(icyl-1)
      CJcylg(icyl)=CJcylg(icyl-1)
      Rcylm=Rcyld(icyl-1)
      endif
      sigma=0.d0
      do jcyl=0,NzDAT
      Zcyl=Zd(jcyl)
C- here dp
C
C               if (jcyl==1) then;
C                  write(*,*)' Zcyl=',Zcyl;
C                  call mypause;
C               endif;
C*
      phia=phiRcylZ(icyl,jcyl)
C- here phia on cylindrical grid
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      rhoRcylZ(jcyl,icyl)=rho
C- 1st index for z
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-             call mypause;
      dRcyl=Rcyl-Rcylm
      if(jcyl.LT.NzDAT)then
      dV=2.d0*((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- multiplied by 2 to account for 'southern' hemisphere
      sigma=sigma+2.d0*rho*dZcyl
C- multiplied by 2 to account for 'southern' hemisphere
      else
      dV=((Rcyl**2-Rcylm**2)+dRcyl**2)*dZcyl*pi
C- not multiplied by 2 since dZcyl is on border
      sigma=sigma+rho*dZcyl
C- not multiplied by 2
      endif
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
C-  	       rhodV(jcyl)=2.d0*pi*rho*(Rcyl**2-Rcylm**2);
      	CM=CM+rho*dV
      CMcylg(icyl)=CMcylg(icyl)+rho*dV
      CJcylg(icyl)=CJcylg(icyl)+rho*vphi*Rcyla*dV
      	CJ=CJ+rho*vphi*Rcyla*dV
      	CI=CI+rho*Rcyla**2*dV
C-                P=EOS(phiu*(phiRcylZ(icyl,jcyl)),4)/(rho0*phiu);
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
C-  	       W=W+rho*phiRcylZ(icyl,jcyl)/2d0*dV;
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      enddo
C- loop on z
      Call GillMiller(Zd,rhoRcylZ(0,icyl),NzDAT+1,ANS,ERROR,IFAIL)
      sigmaGM=2.d0*ANS
C- multiplied by 2 to account for 'southern' hemisphere
      sigma2pir(icyl)=2.d0*pi*sigmaGM*Rcyl
      write(8,'(1p,5g20.8,i5,g20.8)')Rcyl,CMcylg(icyl),sigma,sigmaGM,
     *	ERROR,IFAIL,CJcylg(icyl)
      enddo
C- loop on Rcyl
      Call GillMiller(Rcyld,sigma2pir,NrDAT+1,CMGM,ERROR,IFAIL)
      write(8,'(a,1p,3g20.8,i5)')'# CM CMGM ERROR IFAIL=',CM,CMGM,ERROR,
     *	IFAIL
      do icyl=0,NrDAT
      mcyl(icyl)=CMcylg(icyl)/CMcylg(NrDAT)
      enddo
      VT=abs(2d0*T/W+one+3d0*Pip/W)
      write(*,*)' output controls'
      write(8,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *	,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one
     *,VT
      write(*,2)'cntr',rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,V*Re**3
     *	,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI,T/abs(W),Pip/Eth+one
     *,VT
      if(Lpolytrope)write(8,2)'poly',rho0,r(IB)**delta/alphaPolytrope,
     *	Re/alphaPolytrope,CM*(Re/alphaPolytrope)**3/(4.d0*pi),V*(Re/al
     *phaPolytrope)**3,CJ*sqrt(G_N)*Re**5*rho0**1.5d0/1d50,CI*Re**5
     */(alphaPolytrope**5),
C- *rho0),
     *T/abs(W),Pip/Eth+one,VT
C-        pause;
      LinitMcyl=.false.
C- needed for jmscf -- 1st model is rigid, 2nd jmscf
C------- '<--Leaving  Node %_OutputResults_controls_rhoOnCylGrid:'
C------- '<--Leaving  Node %_OutputResults_controls:'
C------- '<--Leaving  Node %_OutputResults:'
      enddo
C------- '<--Leaving  Node %_tableIB:'
      else
C:  -- IB loop decreasing polar radius *
C------- '-->Entering Node %_loopIB:'
      IB=IA
      do while((IB.GT.3))
C  .and. rotlaw=='alpha') .or.
C                    (IB>0.6*IA .and. rotlaw=='rigid'));*
C- assuming IB<=9999
      IF(.NOT.(IB.LT.10))GOTO 09987
      write(str,'(a,I1)')'000',IB
      GOTO 09984
09987 CONTINUE
      IF(.NOT.(IB.LT.100))GOTO 09986
      write(str,'(a,I2)')'00',IB
      GOTO 09984
09986 CONTINUE
      IF(.NOT.(IB.LT.1000))GOTO 09985
      write(str,'(a,I3)')'0',IB
      GOTO 09984
09985 CONTINUE
      write(str,'(I4)')IB
09984 CONTINUE
      open(26+IB,file='phi'//str(1:len_trim(str))//'.res')
C:  iterate phi etc. for lower IB, which means
C                             a flatter configuration, so faster rotation *
C------- '-->Entering Node %_FindPhi:'
      write(*,*)' entering _FindPhi'
C-            call mypause;
      LogInit=.true.
      write(*,'(a,1p,g15.5,i7)')'rho0, IB',rho0,IB
      write(*,'(a,1p,3g15.5)')'phi(IA,JMax),phi(I0,JMax) =',phi(IA,JMax)
     *	,phi(I0,JMax)
C-           'phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)))=',phi(IA,JMax),phi(I0,JMax),Fpsi(one,rotlaw(1:len_trim(rotlaw)));
C- write of Fpsi is forbidden if there is write in it
      write(*,'(2a)')' rotlaw=',rotlaw(1:len_trim(rotlaw))
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)
     *	*EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      write(*,'(a,1p,3g15.5)')'phiu,EOS(rho0,3),C=',phiu,EOS(rho0,3),C
C-            call mypause;
      n_max=300
!      epsMain=1.d-12
!      epsMain=1.d-9
	 epsMain=1.d-8

C- to data file
      LogIter=.true.
      nStep=1
      do while(LogIter.and.(nStep.LE.n_max))
      C0=C
      Callnew_phi(
C-r,alpha,theta,dtheta, -- now in module cylinderGrid
     *C,LogInit,LogRapid,rotlaw(1:len_trim(rotlaw)))
C- last argument is the rotation law
C              if(C>C0)then;
C                C=C0+min(abs(C-C0),0.03d0*abs(C0));
C              else;
C                C=C0-min(abs(C-C0),0.03d0*abs(C0));
C              endif;*
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18
      if(Dabs((C0-C)/C).LT.epsMain)LogIter=.false.
      write(*,*)' nStep,C,phi(I0,JMax):',nStep,C,phi(I0,JMax)
      nStep=nStep+1
      enddo
C------- '<--Leaving  Node %_FindPhi:'
      if(IB.le.0.2*IA.or.rotlaw.EQ.'rigid')then
      IB=IB-nint(0.01*IA)
      else
      if(IB.le.0.4*IA)then
      IB=IB-2*nint(0.01*IA)
      else
      if(IB.le.0.6*IA)then
      IB=IB-4*nint(0.01*IA)
      else
      IB=IB-8*nint(0.01*IA)
      endif
      endif
      endif
C        do i0=0,IMax0;
C                 rcyl=float(i0)/IMax0;
C                 im=powrsmart(rcyl,one/delta)/dr;
C                 ip=im+1;
C                 rcylm=im*dr;
C                 rcylp=ip*dr;
C                 xm=(rcylp-powrsmart(rcyl,one/delta))/dr;
C                 xp=(powrsmart(rcyl,one/delta)-rcylm)/dr;
C                 phiOut=xm*phi(im,JMax)+xp*phi(ip,JMax);
C                 write(26+IB,'(1p,3e15.4)')xm, xp, phiOut;
C               enddo;
C             *
      enddo
C------- '<--Leaving  Node %_loopIB:'
      endif
C-**      open(16,file='phi');
C-**      write(16,*) phi0;
C-**      close(16);

       Return
C !!!  23.05.2018     stop
C:  *
1     Format(i5,4(1pe10.2),1pe28.20)
2     Format(a,1p,2g17.8,e17.8,g17.8,12g15.6,I5,g15.6)
22    Format(I5,1p,2g17.8,e17.8,g17.8,14g15.6,I5,g15.6)
C-
99999 FORMAT(1X,8F9.4)
99998 FORMAT(F5.2)
99997 FORMAT(1X,I16,F12.4,I11,F12.4)
99996 FORMAT(1X,I39,F12.4)
89999 FORMAT(1X,A,I4)
89998 FORMAT(1X,I4,1p,g20.3)
89997 FORMAT(5X,I4,1p,g12.1,g17.6)
89996 FORMAT(1X,F4.1,A,1p,g12.1,g17.6)
89995 FORMAT(1X,A,I2,A)
      end
C:  SUBROUTINE GillMiller(X,Y,N,ANS,ER,IFAIL)
C              for numerical integration of table data *
      SUBROUTINE GillMiller(X,Y,N,ANS,ERROR,IFAIL)
      DOUBLEPRECISION ANS,ERROR
      INTEGER IFAIL,N
C-     .. Array Arguments ..
      DOUBLEPRECISION X(N),Y(N)
      IFAIL=1
      CALL D01GAF(X,Y,N,ANS,ERROR,IFAIL)
      GOTO(09983,09982,09981,09980),IFAIL+1
      stop' D01GAF produced wrong IFAIL'
      GOTO 09979
09983 CONTINUE
      if(abs(ANS).GT.0.d0)then
      if(abs(ERROR/ANS).GT.0.1d0)WRITE(*,'(1X,A,1p,G12.4,A,G12.4)')'
     *	Integral = ',ANS,'     Estimated error = ',ERROR
      endif
      GOTO 09979
09982 CONTINUE
      WRITE(*,*)'Less than 4 points supplied'
      GOTO 09979
09981 CONTINUE
      WRITE(*,*)'Points not in increasing or decreasing order'
      GOTO 09979
09980 CONTINUE
      WRITE(*,*)'Points not all distinct'
09979 CONTINUE
      return
      end
C:  *
      SUBROUTINE D01GAF(X,Y,N,ANS,ER,IFAIL)
C-
C-     THIS SUBROUTINE INTEGRATES A FUNCTION (Y) SPECIFIED
C-     NUMERICALLY AT N POINTS (X), WHERE N IS AT LEAST 4,
C-     OVER THE RANGE X(1) TO X(N).  THE POINTS NEED NOT BE
C-     EQUALLY SPACED, BUT SHOULD BE DISTINCT AND IN ASCENDING
C-     OR DESCENDING ORDER.  AN ERROR ESTIMATE IS RETURNED.
C-     THE METHOD IS DUE TO GILL AND MILLER.
C-
C-     NAG COPYRIGHT 1975
C-     .. Scalar Arguments ..
      DOUBLEPRECISION ANS,ER
      INTEGER IFAIL,N
C-     .. Array Arguments ..
      DOUBLEPRECISION X(N),Y(N)
C-     .. Local Scalars ..
      DOUBLEPRECISION C,D1,D2,D3,H1,H2,H3,H4,R1,R2,R3,R4,S
      INTEGER I,NN
C-     .. Executable Statements ..
      ANS=0.0D0
      ER=0.0D0
      IF(N.GE.4)GOTO20
      IFAIL=1
      RETURN
C-
C-     CHECK POINTS ARE STRICTLY INCREASING OR DECREASING
C-
20    H2=X(2)-X(1)
      DO80I=3,N
      H3=X(I)-X(I-1)
      IF(H2*H3)40,60,80
40    IFAIL=2
      RETURN
60    IFAIL=3
      RETURN
80    CONTINUE
C-
C-     INTEGRATE OVER INITIAL INTERVAL
C-
      D3=(Y(2)-Y(1))/H2
      H3=X(3)-X(2)
      D1=(Y(3)-Y(2))/H3
      H1=H2+H3
      D2=(D1-D3)/H1
      H4=X(4)-X(3)
      R1=(Y(4)-Y(3))/H4
      R2=(R1-D1)/(H4+H3)
      H1=H1+H4
      R3=(R2-D2)/H1
      ANS=H2*(Y(1)+H2*(D3/2.0D0-H2*(D2/6.0D0-(H2+2.0D0*H3)*R3/12.0D0)))
      S=-(H2**3)*(H2*(3.0D0*H2+5.0D0*H4)+10.0D0*H3*H1)/60.0D0
      R4=0.0D0
C-       write(*,'(a,1pe12.3)')' init ans:',ans;
C-
C-     INTEGRATE OVER CENTRAL PORTION OF RANGE
C-
      NN=N-1
      DO120I=3,NN
      ANS=ANS+H3*((Y(I)+Y(I-1))/2.0D0-H3*H3*(D2+R2+(H2-H4)*R3)/12.0D0)
C-          write(*,'(a,i5,1pe12.3)')' i ans:',i,ans;
      C=H3**3*(2.0D0*H3*H3+5.0D0*(H3*(H4+H2)+2.0D0*H4*H2))/120.0D0
      ER=ER+(C+S)*R4
      IF(I.NE.3)S=C
      IF(I.EQ.3)S=S+2.0D0*C
      IF(I-N+1)100,140,100
100   H1=H2
      H2=H3
      H3=H4
      D1=R1
      D2=R2
      D3=R3
      H4=X(I+2)-X(I+1)
      R1=(Y(I+2)-Y(I+1))/H4
      R4=H4+H3
      R2=(R1-D1)/R4
      R4=R4+H2
      R3=(R2-D2)/R4
      R4=R4+H1
      R4=(R3-D3)/R4
120   CONTINUE
C-
C-     INTEGRATE OVER FINAL INTERVAL
C-
140   CONTINUE
      ANS=ANS+H4*(Y(N)-H4*(R1/2.0D0+H4*(R2/6.0D0+(2.0D0*H3+H4)*
     *	R3/12.0D0)))
C-       write(*,'(a,1pe12.3)')' final ans:',ans;
      ER=ER-H4**3*R4*(H4*(3.0D0*H4+5.0D0*H2)+10.0D0*H3*(H2+H3+H4))/
     *	60.0D0+S*R4
      ANS=ANS+ER
      IFAIL=0
      RETURN
      END
